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FOREWORD 


This interim report summarizes the efforts and accomplishments of 
the subject contract during the period 28 June 1974 through June 1975. The 
study was conducted for the NASA-Dewis Research Center, Cleveland, Ohio, 
by personnel in the Computational Mechanics Section, Lockheed-Huntsville 
Research h Engineering Center, Huntsville, Alabama. The NASA-LeRC 
Project Engineer is Dr. C.C* Chamis, Mail Stop 49-3, 

S. T. K. Chan and C. H. Lee were the principal investigators for the 
study. During the early stage of this investigation, J. N. Reddy was also 
involved and contributed to the mechanics of high velocity impact. The study 
was supervised initially by M. R. Brashears and later by B. H. Shirley. 

Work was begun on this contract on 28 June 1974 and all technical work 
is to be completed by April 1976. Efforts on this contract are being directed 
to the numerical solution of the three-dimensional high velocity impact problem 
based on the hydroelasto- viscoplastic formulation. Provisions will also be 
made to hook up the developed program with existing elements in the NASTRAN 
program through a coupling of Eulerian and Lagrangian modes. 
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AB'STRACT 

A finite element algorithm for^ solving unsteady, three-dimensional 
high velocity impact problems is presented. The computer program is based 
on the^Eulerian hydroelastp-viscoprastic formulatidn and the utilijsatlon of the 

A 

theorem of weak solutions so that the entropy condition is satisfied auto- 
matically. The equations to be solved consist qf^ conservation of mass, 

B i, '! ' 

momentum, and energy, equation of State, and appropriate constitutive equa- 
tions. The solution technique is a time -dependent iinite element analysis 
utilizing three-dimensional isoparametric elements, in conjunction with a 
generalized two-step time integration scheme. Thl^ developed code is demon- 
strated by solving onfe -dimensional as well as three-dimensional impact prob- 
lems for both the inviscid hydrodynamic model and the hydroelastto -viscoplastic 
model. j 
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1. INTRODUCTION 


... ■: ... ■ 

Over the last two decades considerable interest has been shown iiri the 

study of high velocity impact problems. This is primarily due to the need^" 
i ® ; 

for the developmei^t of faster projectiles and the ‘information concerning the 

meteoroid hazard to these projectUes. Current , .interest in high velocity im- 
pact studies is larg^ely due to the concern over the impact of interplanetary 
debris on space vehicles. However, these research results have other 
industrial applications; for instance, in the high t^elocity impact of objects 
on turbine blades, ile., rocks, birds or metals. ' 

I 

Motivated by military applications, seven symposiums were held oj;^ 
the subject during 1955-65 [1-7] under the sponsorship of the Army, Air Force 
and Navy. The bulj: of the material presented aCthese symposiums con s^iod 

of experimental results and very little appeared'on the theories explaining the 

% 

high velocity impact phenomena. During the late 1960s efforts were 
toward formulating ® re alls tic theoiles to explain the complex dynamic, an^ 
mechanical response of materials^in hypervelocity impact [8-12] . 


j The term ” hype r velocity" (or high velocity) refers to the impact vq;^#- 
dity jtegime in which the maximum stress developed by the primary s hook 
wav^^greatly exceeds the material strengths of b6th the target and., the 
jectil^e. During the early stages of the impact process, the target aiid the ^ 
projC,ctile behave essentially as compressible fluids and, c onsequently u 
.seye^al researchers [13-16] have employed pure hydrodynamic mdde^ls to ^M 
analyjse the hypervelocity impact problems. However, these shock vstresse-)^ 
dec^;ig very rapidly as the wave propagates away from the impact point dnd ■* 
reacl( values comparable to the material yield strengths. Fi’om this point ^ 
onwai^d the material strengths become important in determining the stress^ 
and^^elocity fields in the target material and, consequently, due to more 
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involved material response to these rapidly applied stresses, the process 
gets more complicated. Hence, the purely hydrodynamic model cannot ade- 
quately describe the complex physical phenomena [17] and a more realistic 
model must be employed to account for all aspects of the high velocity impact 
phenomena. Several researchers have included the strength eff<^cts (e.g., 
see [18]) and there appeared many computer codes using more realistic 
hydrodynamic -elastic -vis coplastic models, [lO] and [19-23], The numerical 
techniques were largely based on the finite difference method, and they appear 
to be suited for relatively simple (and in particular to hydrodynamic) problems. 

c 

Since a hydrodynamic - elastic -vis coplastic model is basically a structural 
model, it would be appropriate to use structural techniques to solve the problem. 
Also, during the impact process, the geometry involved is generally very com- 
plicated, which would necessarily require a versatile and flexible technique in 
order to treat it accurately and realistically. The finite element method, which 
has proved to be highly successful in analyzing structural problems, is considered 
such a candidate. Leimbach and Prozan [24] have shown the superiority of the 
finite element method over the finite difference method for a simple impact 
problem, although their model is rather simple and unable to predict the actual 
dynamic a.1 response of materials in the impact region. 

The general objective of the present study is to develop a three-dimensional 
finite element code to analyze structural components subjected to high velocity 
impact, using governing equations based on Eulerian and Lagrangian formulations. 
In particular, the impact point with its vicinity is to be represented by an Eulerian 
hydrodynamic -elasto-viscoplastic model, while the remaining structural com- 
ponents are to be analyzed with existing computer programs for structural anal- 
ysis such as NASTRAN. To bridge the gap between the Eulerian mode used for 
the impact region and the Lagrangian mode generally employed in structural 
analysis, a finite element code based on coupled Eulerian and Lagrangian form- 
ulation is being developed to match the two solutions. The computer program 
will, at the end, be capable of handling static and dynamic response of large 
deformations, anisotropic material behavior, plastic yielding and material 
fracture. 
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---^^This report contains mainly 1;he development of a finite element Gx>m-i 

. • . ■■ . • 

pater program for the numerical solution of three-dimensional high; veTocify 

impajft problems, based on the Eulerian hydroelasto -viscoplastic formdla'tixbn . 


he. equations to be solved consist of conservation of mass, momentuni, arid 
AY f equation of state, and appropriate constitutive equations. The prirfi- 


i.fci^nvaHables . i,ei^ density, momentum, total en-ergy, and stresses, are used 


, . as p^mary unknowns in the computations. The solution technique deyeloped 
in tb^ study is a time -dependent finite element analysis utilizing tbfe^e^^ 
dimensional isoparametric elemeiits. The finite' element analog of the govern- 
ing equations is constructed as a Consequence of the theorem of weak solutions, 
SO that the entropy condition can be satisfied automatically in the formulation. 

As an integrated part of the algorithm, a generalized tv'-o-steo, time-splitting 

® V ' • ' 

finite element scheme is proposed^ to remedy the^jpumerical instabilities during 

time marching. , • 


Following this introduction, a description of ^the Eulerian, Lagrangian, 
and coupled Eulerian -Lagrangian formulations are presented in Section 2. 

The latter formulation is believed to be useful in the latter stages of develop- 
mentj In Section 3^ the mechanics of high velocity impact are discussed and 
the_^verning equations are presented in the conservative (divergence) form, 
alon^with the specific form of the, equation of state for metals. A. general^ 
dis.ci^sion of initial and boundary conditions is also presented. Section 4 
Heyofled to the finite element formulations investigated earlier, utilizing epn- 
venti2>nal methods of weighted residuals, for governing equations with. Eul.e-r ian 
desc^’ription. In addition to problem formulation, element desc riptlonx rand 
^I^UIrl^ical integrations are also discussed. Two time integration ©.cLemesl 
nam^y, the implicit finite difference and the Galerkin in time, arc dcsc;ribc.*d 
.at.th^ end of this section. Presented in Section 5 is an improved finite ele- 
algorithm based on the Eulerian hydroelasto-viscoplastic formulation, 
nsr&K the theorem of weak solutions, together with a generalized two-step 
time i splitting scheme. Free surface considerations and large system-equa- 
tion Solver are also discussed and will be integrated into the final version of 
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the computer code. Section 6 summarizes numerical results and findings on 
a number of test problems. They include the heat conduction problem in solids, 
impact problems computed with an in viscid hydrodynamic model, and a problem 
computed with the hydroelasto -viscoplastic model. Summarized in Section 7 
are the findings up-to-date together with further investigations to be conducted 
during the remaining period. 
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2. VARIOUS DESCRIPTIONS AND FORMULATIONS 


Consider an open (i.e., not including its boundary') bounded region f2 
at a time t = 0 in three-dimensional Euclidean space with its boundary 0S7 
(see Fig. 2-1), The7 union of and its boundary BQis the complete region 
(occupied by a body; at time t = 0) and is denoted hy 




^ Q U BQ 
o 


where U denotes tl^e union of twd'sets. Thus, in Fig. 2-1, denotes the urjion 
of and its boundary Note that the region has two boundaries: one 

external boundary and one internal (the interface) boundary 9Q^. 



Fig. 2-1 - A Three-Dimensional Region at t - 0 
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Suppose that we are given, at time t = 0, the state of the fluid (or solid) 

occupying the region the external forces acting on the boundary, and the 

boundary conditions. We are then required to determine the state of the fluid 

and the shape of the -region at a subsequent time, t = T. Since we plan to 

analyze the impact point and the surrounding region, j, separated from the 

remaining region, it is convenient to consider two different fluids (or materials) 

occupying two different portions, 17, and 17^, of 17 ; that is 

1^0 

u u 

S- 

The surface 817^ is the material interface of the two regions, and it moves, 
as time advances, in such a manner that the pressure (or stresses) and velo- 
city components are continuous. 

2.1 EULER AND LAGRANGE DESCRIPTIONS 

There are two basic descriptions well known in continuum mechanics 
with respect to which the governing equations can be derived and computations 
can be carried out. In the Lagrange description, mostly used in solid and 
structural problems, a coordinate system is fixed in the body or configuration 
17^ to be studied. The deformation of the projectile -target configuration is 
then measured with respect to this deformed configuration. Consequently, 
the positions of the boundary and of the material interface 912^ are auto- 
matically determined. The description also permits the use of constitutive 
relations for the material in which the stress history of each portion of the 
body is taken into account. However, the Lagrange description is totally 
unsatisfactory for calculating a flow in which turbulence develops or in which 
new material interfaces develop. In this case the nodal points of the mesh 
will attempt to follow the motion and the material particles which were initially 
adjacent to each other in no longer remain so as the turbulence develops. 

In Eulerian description, mostly used in fluid mechanic problems, the 
coordinate system is fixed in the space rather than in the body or configura- 
tion, and the calculations follow the material point that happens to be in a 
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given location at that particular time. In this case the large distortions do 
not cause any problems, however, more than one material cannot be treated 
accurately. The curves approximating and DH ^ move with the body and 
therefore create irregular, time -dependent boundary zones In the fixed 
Eulerian mesh. 

2,2 COUPLED EULERIAN- LAGRANGIAN DESCRIPTION 

' Obviously, Aeither Lagrange description nor Eulerian description alone , 
is ideally suited for the analysis of impact problems. It is both natural 
and desirable to combine the use.pf Eulerian and Lagranglan modes depcndiiTg 
on whether the material is in a ](:luid state (Eulerian mode) or solid state 
(La^randian mode). In doing so, there will be a great deal of flexibility in 

3 - 

approximating the -problem, thus enhancing the solution process regarding 

accuracy and computational efficiency. The idea is very similar to the well 

known sub structuring technique, but now with different descriptions (or modes) 

used in various regions. More specifically, in analyzing the impact problem 

of structural components, the existing NASTRAN program (in Lagrangian 

mode) can be used to advantage for the structure part, while an Eulerian 

*16 

description is necessary for the impact point aifd its vicinity. In the following, 

i. 

we discuss how to couple an Eulerian mode with a Lagrangian mode, in partic- 
ular, the NASTRAN program. - 

Conceptually, the coupling of Eulerian and Lagrangian modes can be ac- 
complished through the use of a Coupled Eulerian-Lagrangian Code (CEL), both 
involving velocities, v, as unknowns [25] . Furthermore, if it is desirable to 
use also the NASTRAN program with displacement, u, as unknowns, another 
coupling of the Lagrange mode with the NASTRA^ program (CLN) must also be 
considered. Figure 2-2 shows a configuration consisting of these different 
descriptions. In the figure, the dotted lines represent the interfaces between 
different modes. The interfaces belong to the Lagrangian mode and are to be 
adjusted as time advances. Therefore, the coupling of an Eulerian mode with 
the NASTRAN mode is to be accomplished through two steps, i.e., the CEL and 
CLN codes. 



2-3 


LMSC-HREC TR D390900 



Fig. 2-2 - A Typical Configuration to be Analyzed Using 
Various Descriptions 


The basic idea of the CEL formulation is to approximate a configuration 
by the combination of Eulerian and Lagrangian subregions, with the boundary 
and interfaces described by Lagrangian lines. The Eulerian mode will con- 
sequently have its boundary prescribed by the Lagrangian calculations. Thus 
the Eulerian calculation reduces to one based on a fixed mesh but having a 
prescribed moving boundary to reflect the coupling effects. 
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Consequently, the calculations that are made at ea^h time step will involve 

three parts: JLagrangian calculations, Eulerian calculations, and a calculation 

which couples the Eulerian and Lagrangian modes by keeping track of the 

moving interfaces and the matching of pressure, velocities, etc., along this 

interface. Suppose that we know the solution of field variables (energy, density, 

velocities, etc.) at the n^^ time step, t^, and wish to compute the solution at 

n+1 

the next time step, t^ . Since we know the position of the Lagrangian mesh, 

L^, arid the Euleriah'mesh, E*', with its boundary subject to move, the calcu- 
lations for the next time step can proceed in the following way. 

.1 ■<? : 

The first calculation uses the known (t=t^) state of the Lagrangian region 
and the pressure acting on its boundary to solve th^biequations in Lagrangian 
description. The soflition gives us the t = t^^^ state in this region and also a 
new mesh position, V 

The boundary of the Eulerian region is then updated, using results from 
L^^\ and computation is performed to seek for the t”^^ state of the fluid in 
the region E. This iS considered as the Eulerian phase of the CEL calculation. 

n+1 

Having determined the t state of the Eulerian region, the second phase 
of coupling is to determine the pressure to act on the boundary of the Lagrange 
region L^^^. We have thus advanced all of the field variables and grid positions 
to their values at the next time step, which completes one basic cycle of the CEL 
calculations. 

Similar logic and procedures can be followed for the coupling of L (v) with 
NASTRAN. Herein, the conversion between velocities and displacements can be 
done through integration in time . Another consideration is the generation or 
suppression of nodal unknowns on the interface, which can be accomplished using * 
the approximating polynomials. 

Another possible way is to couple the Eulerian mode with NASTRAN directly 
(CEN) as shown in Fig. 2-3 on the following page. In this case, no formulation for 
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the L(v) region is necessary, but the coupling must be done through the E-N 
interface. This is probably somewhat more complicated because the coupling 
involves simultaneously two kinds of descriptions in terms of different field 
variables. However, in light of the possible savings in computational aspect, 
this approach is also worth investigating. 



Fig. Z-3 ~ Coupling of Eulerian Mode with NASTRAN 


2.3 LOCAL REPRESENTATION OF DYNAMICAL SYSTEMS UNDER AN 
ARBITRARY FRAME 

In order to construct a code coupling the Eulerian-Lagrangian modes, 
it would be more convenient if the field equations are presented in their local 
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form under an arbitrary frame of reference. More specifically, lot x = 

be a Cartesian coordinate system moving with a certain velocity rela- 
tive to the motion of the physical field, and let X = (Xj^, X^, X^) be the refer- 
ence system. Then, the two coordinate systems can be related as 


X = X (X, t) 


(2.1J 


When we know such a relation for every material particle in the medium, we 
say that the history of the motion is known at time t. We assume that the 
motion is continuous, single valuedtand that (2,1) can be inverted to give the 
initial; position or material coordinates, X. A necessary and sufficient condi- 
tion for the inverse to exist is that the Jacobian should not vanish: 


J = 


8x. 


8X. 

J 


;^0, 0<J<oc vt>0 


( 2 . 2 ) 


Let Q denotes the velocity of material element relative to the x system, 


i.e., 


Cl ^ 


d X 

"dT 


= V « X 


(2.3) 


where ^ is tjie velocity of the material element, and x = 8x/0t is the local time 
rate of change of the x system. If F represents a certain physical quantity 
satisfying a dynamical system, then^the total time derivative follows the equality 


^ - 
dt 


8 t 


+ (f2 • V) F 


(2.4) 


Clearly, when x = 0, the systems x-and X are idential, and x represents the 
Eulerian system; when, x « v , on the other hand, it represents the Lagrangian 
systern. 
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Following the implicit differentiation rule, the divergence of the local 
rate of change of x system can be related to the Jacobian as 

V-x=i/j (2.5) 

thus, (2.3) yields 

V-S=V.v-j/j (2.6) 


Suppose now that there is a dynamic system with physical quantity F 
satisfying the equation 


dt 


G (x, t) 


(2.7) 


In particular, the continuity equation in x system is readily obtained as 


^ + pv-v = 0 (2.8) 

at 'v.' 

Then, the dynamic system (2.7) can be deduced by applying (2.3), (2.6) and (2.8) 
into the following form: 

+ p fa + V *(QpF) = pG (2.9) 

r) t ^ 


where 

A = j/j 

It. is clear that the dynamic system (2.9) would reduce to the Eulerian form 
when X = 0; and to the Eagrangian form when 2=0. 
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3. MECHANICS OF HIGH VELOCITY IMPACT 


The high velocity impact of a projectile with a solid target results in 
an extremely complex phenomenon. A complete description of this problem 

, V-. 

would involve considerations of all phases in the theory of continuum me- 
chanics. This includes not only the compressible fluid flow, dynamics of 
elasticity and plasticity, but also other behaviors such as melting and solidi- 
fication, vaporizatibh and condensation, and the kin'etics of phase change. 

Several models have been proposed for the various stages of the high 
velocity impact problems (see, e.g., [10] , [ 13-16] j [ 19-23] ). Two models 
corresponding to different stages, namely, the inviscid hydrodynamic model 
and a hydroelasto-viscoplastic model, are proposed and analyzed in this study. 
The following discusses the mechanics of high velocity impact. 


: . 

3.1 IMPACT PROCESS 

\ 


"The analysis of high velocity .projectile mechanics can be divided 
into two parts: (1) dynamic behavi(j)r of the projectile and target during the 
penetration process^ and (Z) structural response df the target after the pene- 
tration process is completed. g, ^ 

ft 

& 


During the short period of time in which the projectile contacts the 
target, a plane shock wave is generated in the projectile as well. as in the 
target. The pressure behind these shock fronts is the largest pressure that 
exists throughout the entire impact process. Due to the high pressure the 
material strength of the target can be ignored, and the material can be 
assumed to behave essentially as an inviscid, compressible fluid^' The shock 
waves generated in the projectile and target travel away from the interface 
(see Fig. 3-1). If the projectile is finite (in diameter), rarefaction waves will 
be ge|ierated and transmitted toward the axis of symmetry. 



It- 


j 





Fig. 3-1 - High Velocity Impact Process 
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Formulation of rarefaction waves results in the ejection of target and 
projectile material particles. Moreover, the rarefaction waves weaken the 
shock waves and change its plane shape to approximately spherical. The 
strength of the shock continues to decrease due to the spherical attenuation 
and additional rarefaction waves, and the influence of the material strength 
and the strain rate effects must be taken into account. 

■ ^ 

3.2 RANKINE-HUGGNIOT RELAT^lONS 

As aforementioned, the material can be considered as an inviscid, com- 
pressible fluid in the short period right after the impact. The problem, in- 
cluding the geometric development of the impact shock, can then be treated 
as an unsteady, supersonic flow resiembiing a moving shock. 


Assuming that the hemispherical shockwave profile is steady in time, 

the Rankine-Hugoniot relations relating the pressure, P, internal energy, r , 

and the density, p, behind the shock to the same quantities in frorit of the shod. 

are applied. These ‘equations express the conservation of mass, momentum, 

and energy in termsfcof the shock velocity, v , and particle velocity, v : 

~ s p 


P v 
*^o s 


= p(v - v^) & 

s p 


(3.1) 


P - P 


- p v v 
*^o s p 


(3.2) 


2n 




p V 

^o s 


= P V. 
o s 


(3.3) 


where the subscript o refers^to the initial (or undisturbed) values. The 

product p V is called the shock impedance . Using (3.1) and (3.2), (3.3) 
o s 

can be written as 


e - 


= 


+ p )(^ -^) 
° Vo p> 


(3.4) 
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Solving for v and v from (3.1) and (3.2), we obtain 

B p 

• "P ■ 'P/P„-1I (3.51 

The initial pressure, P^, is generally very small and can be neglected in 
cbmparison with P. If shock velocity and particle velocity are known at a 
point on the shock front, the pressure can be computed using (3.2). One 
more relation is needed to solve the flow across a shock, namely, the equa- 
tion of state; a specific form of the equation of state is given in Subsection 
3.4. 


3.3 BASIC EQUATIONS 

As pointed out previously, an inviscid hydrodynamic model is a good 
approximation in the early stages of the high velocity impact when the pres- 
sures developed are much larger than the shear strength of the material. 

It ceases to describe the phenomenon adequately in the later stages of the 
impact when the influence of the shearing strength cannot be neglected. 
Therefore, the elastic, viscous and plastic effects of the material must be 
included. The equations governing the dynamic behavior of elasto- 
viscoplastic materials are described below. 

3,3.1 Conservation Equations 

The conservation of mass, momentum and energy leads to the following 
equations by letting A = 0 and 2 “ X 

It + ^ 

1 

-;^(pv. ) + (pv- V. ) = 

3t l ' Bx. j 1 dx. ' ij 

3 J 


+ pfi (3.7) 
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J 


-ir. + Ph '"i + p® - 8^ 

j 1 


(3.8) 


where p is the density, is the velocity in the x. -direction, is the stress 
tensor, is the body force per unit mass, e is the specific total energy 
defined as ^ n 


e = e + ± V. V. 

2 11 


(3.9) 


e being the internal energy per unit mass, S is the rate of internal 
specific heat generation per unit mass, and is the heat flux or the rate of 
heat flow per unit a^ea across the^,surface in thesdirection of its unit outward 
normal, n^. f\.\ 

3.3.2 Hydrodynamic-Elasto- Viscoplastic Constitutive Relations 

■ When the medium under consideration is inviscid, there are no shearing 
stresses. The stress tensor then .becomes 


a. . ^ - P 5. . (3. 1C 

ij U 

*s-. 

where P is the hydrostatic pressure that is independent of orientation, and 
6. . is the Kronecjce^ delta. In general, the stress tensor and the correspond- 
ing ^rain rate ten^r are relatedOto their respective deviatoric -tensor s by 

the following forms^ 

■ (* c- 

S. . = CT. . - i a, , 6. . (3.11 

1) ij 3 kk ij ' 

d.' . = d. . - i d, , 5. . (3.12 

ij ij 3 kk ij 


where S. . is the deviatoric stress tensor, d. . is the strain rate tensor defined 

ij n 

"" /3v. 9v.^ 

d,, = T +^) (3.13) 

J 


ij " 2 \8x 

I 
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and d/j is the deviatoric tensor of strain rate. The deviatoric stress tensor, 

S. in turn, is related to the strain rate tensor d. . through another constitutive 
equation. For linear elastic materials, the constitutive equation for stress 
and deformation rate is given by 




(3.14) 


where is a fourth-order tensor of material parameters. If the stress 

components are symmetric, is also symmetric, i.e., = ^ijik ~ 

Eji^k* general, there are 21 independent elastic constants. For ortho- 

tropic materials, the number of independent constants reduces to 9, and for 
isotropic materials it reduces to 2. The isotropic constitutive relation is 
given by 


S. . = 2/i d. . + X± , 6. , (3.15) 

ij ij kk ij ' ^ 

where X and /i are the Lame’s (or viscous) constants. 

Let us denote the mean of the principal components of stress by 

Q H- ^ (3.16) 

The quantity Q is called the dynamic pressure. The rate of dynamic pres- 
sure dQ/dt may, in general, be decomposed into a thermodynamical reversil)le 
component dP/dt and a dissipative component djr/dt, i.e., 

f=f^f ( 3 . 17 ) 

In general, the rate at which the local thermodynamic equilibrium is attained 
is much greater than the rate at which a disturbance can be propagated. It is 
then reasonably accurate to assume that the local thermodynamic equilibrium 
exists at each instant. Hence, the reversible rate of pressure, dP/dt, is not 
path-dependent, and its integral P follows the equation of state which can be 
expressed as a function of density p and internal energy c, i.e., 

P = P(P, e) (3.18) 
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The dissipative component drr/dt, on the other hand, is generally a path- 
dependent function related to the bulk viscosity and the rate of the change of 
volume. It characterizes the physical dissipative rate of dilatation. For 
isotropic materials, 

n = - a + I M) (3.19) 

The coefficient of bulk viscosity is defined by 

= (\ + I (i) = (P-Q)/dj^j^ (3.20) 

If the volumetric changes of the materials are elastic, the changes of dynamic 
pressure is reversible which implies that dQ/dt =dP/dt. Then the path- 
independent nature of P yields 

Q = P = P(p, e) (3.21) 

Noting that for incompressible materials (dj^j^ = 0), Eq.(3.21) follows 
immediately. Equation (3.21) is also true for compressible fluids for which 

(X + y )i) = 0 (3.22) 

Condition (3.22), known as Stoke 's hypothesis, is a reasonable assumption for 
flow of monatomic gases; however, it is not valid for polyatomic gases or 
liquids, and distinction must be made between the mean stress Q and the 
thermodynamic pressure P, 

For high velocity impact problems, the thermodynamic pressure is very 
high, and as aforementioned, the dynamic response in the material can thus be 
considered as an is.entropic process. This implies that the material under 
high velocity impact can be assumed to possess the elastic changes in volume, 
and (3.21) follows. 
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Unfortunately, there does not exist a constitutive relation that describes 
all aspects (elastic, viscous and plastic) of mechanical behavior in a single 
expression. Constitutive relations of plasticity and viscoelasticity are essen- 
tially dynanlical in nature. The constitutive relation in classical plasticity 
involves tensors of stress and rate of deformation, and describes rigid, per- 
fectly plastic behavior. When elastic effects are to be considered in the 
analysis, this relation applies to the plastic^ part of the rate of deformation 
tensor. Similarly, the constitutive equation for viscoelastic material involves 
stress, elastic rate of deformation and the rate of deformation of a viscous 
fluid. In both cases, the elastic component of the rate of deformation tensor 
is usually written as a function of stress rate. Although the stress tensor is 
objective (axiom of objectivity requires that if a stressed continuum performs 
a rigid body motion and the stress field is independent of time when referred 
to a coordinate system attached to and moving with the material, the stress 
rate must vanish identically), the stress -rate tensor is not. Therefore, the 
stress rate da. ./dt should not occur in this form in the constitutive relation. 

An objective tensor containing the stress rate must be defined. Several ob- 
jective tensors containing the stress rate can be constructed. One such tensor 
is due to Jaumann (see [ Z6] ) and is given by 




ki 


dt 


+ a. 


km ^mi 


- cr „ w, 
mf km 


(3.23) 


where 



(3.24) 


The tensor a^^^ is called a stress flux . Other stress fluxes may be obtained 
by adding objective tensors such as ^mf right-hand side of (3.23): 


da 


ki 


ki 

dt 


a . V . + 
mi m, k 


cr, V . 
km m, i 



da. 


ki 

dt 


- a V - o, V 

mi k, m km i,m 


(3.25) 


da 




’ki 


dt ^mi '^k, 


m 


a, V , 
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The stress flux measures the rate of change of the stress components 
with respect to a rectangular Cartesian system that participates in the rota- 
tion of the material, and a.j = 0 implies that the invariants of stress tensor 
are stationary. 


In general, the rate of deformation tensor d^j can be split into two parts, 


i.e., 


d.. = d^ 


(3.26) 


where dfj is the elastic component and d^? is the viscoplastic component of 

the deformation rate tensor. Here it is assumed that the medium is initially 

' , 0 I 

unstressed. For linear elastic material, the elastic part d. . is related to S. . 

by Eq. (3,14), and for linear elastic -viscoplastic rinaterial d^^ is related to 

the stress flux S. . as 
iJ 


“ ^ijkf 


(3.27) 


in the elastic region, . Here the coefficient matrix depends on 

Eq.(3.14). For isotro]f)ic material satisfying constitutive relation (3.15), Eq. 
(3,27) can be deducedHurther as , 


A 

s. . = 
u 


d S. . . , 

— + s. 0) . - s . u). 

dt im mj mj im 


= 2m d. . 

1 j w 


(3.28) 


Plastic behavior is assumed to take place when a certain function of 

I 

deviatoric stresses vanishes. This function is callbd the yield function. The 
yield function is constructed based on the following assumptions; 


The yield surface is convex (or smooth). 

The plastic component of the deformation rate tensor 
is nornnal to the yield Surface at a smooth point. 
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If a yield condition is given by 

F(cr.j.K) = - Y(»f) = 0 (3.29) 

•A 

with F < 0 denoting the purely elastic region, then the deformation rate tensor 
will be a function of positive values of f. Here Y is a yield stress and ^ is a 
history dependent hardening (or softening) parameter. We now introduce the 
notion of ‘’plastic potential“^{a^^) defined as 


V - 

U 


d: 


(3.30) 


where y is a fluidity parameter which may depend on time, invariants of the 
rate of deformation, etc., and f^ denotes a reference value of f that makes the 
expression non-dimensional. To ensure no viscoplastic flow below the yield 
limit we write 


< 




J ° 

U(fAo). 


F < 0 


A 

F > 0 


(3.31) 


A sufficiently general expression for 0(f/f^) is given by a power law 

(h (3.32) 

Various yield criteria and plastic potentials can be introduced dependino; 
on the material under study. For isotropic materials, for example, there arc^ 
two well-known yield criteria: the Tresca yield criterion assumes that the 

yielding occurs when the greatest difference between any pair of the principal 
stresses cr^, ^.nd reaches a specific value, Y. Usually, Y is taken to be 
the yield stress in uniaxial tension. The von Mises yield criterion for plastic 
yield is 


(ct^ - ^ 


or expressing in terms of 


we have 
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<‘^11 - ^ 22 ^^ + <<^22 - ® 33 >^ + <«^33 r ^ ^ <<^12 + <^13 + ^ ( 3 . 33 ) 

In what follows we shall take a more general yield function that is defined in 
terms of stress invariants as 


f(cr..) = f(Ji,J2, J3) 


where J,, and are the stress invariants, 
' 1 z j j 


J , — o*. J— — y Sip. S.., 

1 11’ 2 2 13 ji’ 3 


T S..S., S. 
3 ij jk ki 


Clearly, von Mises yi.eld criterion, Eq. (3.33), is a special case of Eq. (3.29), 

A 

where F is given by ' 

“ - Y = F(S..) (3. 

where Y is the uniaxial yield stress. 

A *5 

For a special cs^se in which 4^ =; F and Y is independent of we have 
from Eq. (3.30) ^ 


9 Y(f/fo)" »F/8<7. . = y(f/f^)" 8f/8S. . 


and for n^ = 1, and f giVen by Eq. (3,35), this reducers to 

After normalizing the yield surface by the von Mises yield stress, i.e., letting 
F = F/Y, a more general form analogous to (3.36) follows from (3.35), namely 


<7 = y0(F)s../^ 


where 


F = f/Y - 1 

▼ I Y 
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and the plastic potential 

( 0 for F < 0 

0(F) = 

(0(S..) for F > 0 

' ij' 

Thus, the stress flux for this case becomes upon substituting (3.37) into 
(3.28), while noting = 0, 

S.. = ZuidU - y<HF) S /JT;) ' (3.38) 

IJ Ij IJ W 

Here, the material parameter y and the plastic potential 0 have to be found 
from experiments. For example, y and 0(F) of 2024-T3 aluminum are readily 
obtained as 

9 

0 = - 1 for 0 < 0^ 

0 +1 

^ ' ^o ^ ^ ^o 

and y =0.15 sec*\ where 0^ = 10^, F^ = 9 log(0^ + 1), and 6 is a new 
material function which depends on the generalized plastic strain Cp as 

0(Cp) = 0.0003 + 0.0214 - 0.0243 6p for gp < 0.4 

= 0.005 for > 0.4 


Here, the generalized plastic strain is 

t 
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3.3.3 Equation of State 


k 


Since the governing equations contain the pressure (implicitly’* in a..), 
a constitutive relation must be used relating the pressure to the density and 
the specific internal energy. This particular constitutive relation is well 
known as the equation of state. For high velocity impact of solid bodies, 
two equations of state are well known. One of these is developed' by Tillotson 
[27] and the other is developed by Osborne and his associates at Los Alamos 
Scientific Laboratory. 

Tillotson’s Equation of State: Tillotson's equation of state (for either 
compression or expansion regions) is given by y. 


P = 7T (e.Eip) + Afi + o,. 


(3.39a) 


for p > p with 0 < e < € » and 
o s 


P = acp + 


beP 


[d + e/e„n^' 


. + Ap e 


-Pd/n-i) 




(3.39b) 


for p < p with e > e . where 
o s 


p = I rj -1, r] = 


■ifii 


7r(6,p) = ePi 


p/p^, p^ = initial density 


e/ 


and a,b,A,B and are parameters (constant) which depend on the material, 
and 6 is the sublimation energy. , Values of these parameters for some 
materials are given in Table 3-1. 
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Table 3-1 

VALUES OF THE PARAMETERS IN TILLOTSON'S EQUATION 


Parameter 

Aluminum 

Iron 

Copper 

Lead 

p^(gm/cm'^) 

2.702 

7.86 

8.9 

11.34 


a 

0.5 

0.5 

0.5 

0.4 


b 

1.63 

1.5 

1.5 

2.4 


A (mb) 

0.752 

1.28 

1.39 

0.466 


B (mb) 

0.65 

1.05 

1.1 

0.0026 


a 

5.0 

5.0 

5.0 

10.0 


(3 

5.0 

5.0 

5.0 

2.0 


e (mb -cm /gm) 
o ^ 

£ (mb-cm /gm) 
s 

0.05 

0.03 

0.095 

0.0244 

0.325 

0.0138 

0.02 

0.0026 



Note that (3.39) is not defined for certain states, for example when p < 
with e < € , and p > p with e > e . Some investigators have used some 
kind of average of the pressures given by (3.39a) and (3.39b) for these states. 
At this writing Tillotson's equation is not used in the calculations. 

Los Alamos Equation of State ; Another equation of state that is widely 
employed in high velocity impact calculations was derived at Los Alamos 
Scientific Laboratory. The equation is given by 


P(e,p) = 


[Am + 6E (B + 6EC)]/(6E + tp^), M > 0 

[m Aj + 6E (B^ + nBj + 6EC)]/(6E P <0 


(3.40) 


where 

6 = p^/p, p = p/Pq - 1 

A = Aj+pA^, B = + p (B ^ + fi 

C = + p C E = p€ 
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and A A^. B C^, C j and -ip^ are constants that depend on the material. 

Values of these parameters for plexiglas, graphite, aluminum, iron, copper 
and lead are given in Table 3-2. . 


Table 3-2 

VALUES OF THE PARAMETERS IN THE LGfS ALAMOS EQUATION 
OF STATE FOR- VARIOUS MATERIALS 


Parameter 

Plexiglas 

Graphite 

Aluminum 

Iron 

Copper 

Lead 

Po 

1.18 

2.25 

2.702 

7.86 

8.9 

11.34 

Al 

0.006199 

0.1608 

1.1867 

7}.78 

4.9578 

1.4844 

-2 

0.015491 

0.1619 

0.763 

31.18 

3.6884 ' 

1.6765 

B 

o 

0.14766 

0.8866 

3.4448 

9.591 

7.4727 

8.7317 

®i 

0.056H 

0.5140 

1.5451 

15.676 

11.519 

0.96473 

^2 

0.50504 

1.4377 

0.9643 

4.634 

5.5251 ' 

2.6695 

c 

o 

0.5575 

0,5398 

0.43382 

0.3984 

0.39493 

0.27732 

<=i 

0,6151 

0.5960 . 

0.54873 

0.5306 

0.52883 

0.43079 


0.100 

0.500 

1.5 

9.00 

3.60 

3.300 


The parameters are fitted for gram -c entimeter -microsecond system of units. 


A comparison of the pressures computed from Tillotson’s equation of 

state and Los Alamos equation of state is made for aluminum (p = 2.702 

3 3 ^ 

gm/cm ) at e = 0.5 mb-cm /gm (see Fig. 3-2). The equation of state which 

closely agrees with^the experimental data for a given material will be used. 


Hugoniot Equation of State : Ais mentioned earlier, we need to add an 
equation of state to the Rankin e- Hugoniot equations (3,1) through (3.3). For 
instance, Eq. (3.40) may be used and with the internal energy e eliminated 
using the Rankine-Hugoniot jump relation (3.4) to obtain 


P = J 


- 

f 

2 

^2(6) ‘J 

' 

’ j 

2 

£3(6) - J 

* 


A 

i 


=2< 


/Fj(6), 6<1 

/Ej(6), 6>1 


(3.41) 
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where 6 = 

jt 

Ej(6) = d'^'d') [(1 - 6) |(C^ - Cj) 6 + Cj'l - 26] 

Ea'S) = - ^ r<®o - Bj) 6 + Bj] + 6 

s 

£3(6) = Aj (1 - 6) , 

Fj = E36, F2 = E^d - B2U - 6)V2 
F3 = E36 + A2 (1 - 6)^ 

and the constants A^,A 2 , • , are the same as those appearing in (3.40). 

3.4 INITIAL AND BOUNDARY CONDITIONS 

To complete the description of the high velocity impact problem, we 
must include the initial and boundary conditions, which, in general, are 
time dependent. 

Initial Conditions ; At time t = 0 values of all the dependent variables 

(p, V, e, must be specified at the nodal points of the mesh. It is not essen- 

tial to specify all these quantities at the same set of nodal points. For 
instance, in the analysis 6f target alone, the velocities, density and specific 
internal energy are specified at one set of points, and pressure is specified 
at a different set of points, 

\ 

Boundary Conditions : Depending on the type of the boundary (e.g,, 

rigid boundary, free surface, interface, plane of symmetry, etc.), there 
are different kinds of boundary conditions in a problem. At a rigid boundary 
the normal component of thd particle velocity must coincide with the normal 
component of the velocity of the rigid boundary. For a fixed (in time) bound- 
ary, the normal component of the particle velocity must be zero at that 
boundary. A plane of symmetry can be interpreted as a fixed boundary. 
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On a free surface, the total stress must vanish. At an interface (and at a 
contact discontinuity) the total stress and the normal component of particle 
velocity must be continuous, and the density, internal energy and the tan- 
gential component of particle velocity may be discontinuous (jumps). Across 
moving shock fronts the Rankine-Hugoniot relations (3.1) through (3.3) must 
be satisfied. 
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4. BASIC FINITE ELEMENT PROCEDURES 

In this section, we summarise the basic fiMte element procedures as 
applied to solve the impact problem as governed by a set of hydrodynamic 
and, constitutive ec^uations. The formulations are based on the conventional 
methods of weighted residuals, i.e., the Galerkih and least squa;res approches. 
First, the governing equations to be solved are presented. Finite element 
analogues for these equations arc then constructed via the Galerkin and least 
squares approaches. In numerical aspects, the isoparametric elements together 
with numerical'integration are described, and finally, two tiliie -marching 
schemes are presented. The procedures discussed here serve as a reference 
for the development of an improved model described in the next section. Also, 

I 

the numerical results so computed provide data for comparison purposes. 
These methods have been found tc| be quite effective for solving problems 

« * i> 

governed by elliptik: and parabolio equations. Hbwever, for hyperbolic system 
of equations, numerical instabilities were encouhtered, indicating conventional 
methods of weighted residuals cannot be directly applied to solye the impact 
problem. For these reasons, an^^mproved scheme based on the theorem of 
weak solution was "developed, whi6h is described in the next section. 

I 1^ 

& 

4.1 GOVERNING 4:QUATI0NS 

T tl 

With the internal heat generation and heat Tlux set equal to zero, the set 
of equations conserving mass, momentum, and energy becomes 


^(P) + (pv.) = 0 


(4.1) 


J 




(4.2) 
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= -WT. '^i 


(4.3) 


For the purpose of numerical computations, we denote 


pv. = V-, pe = E, pf. -=5. F. 


3v 


3v, 0 V., Ov. 


^ UVk c. vj c. 

1 " ft* _ _ 


8x 


0Xj * 8x^ 


(4.4) 


and rewrite (4.1) through (4.3) in an alternate form 


|e t v,-^(pit 0 p . 0 


(4.5) 


JL 

at 




(4.6) 


W (E) + V. ^ (E) + 0 E = ^ (Oj. -)+ F.v. 


(4.7) 


The stress tensor, a.., can be expressed as the sum of a deviatoric 

^ J 

stress and a dynamic pressure, namely. 


o. . 

13 


S.. + 

13 


I 

3 


a 6 . . 
kk 1} 


(4.8) 


For high velocity impact problems, the first invariant of stress tensor is ex- 


pressed by 


o, 


kk 


- 3P 
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Th’e deviatpTic stress is, in turn, governed by a constitutive equation in the 


form 


as.. as.. e« 

= ^ijkf ■ ®im “mj ^mj “im 


(4.9a) 


For isotropic materials satisfying von Mises criterion, the above equa- 


tion reduces to 




8s. . as.. 


[d.'. - v<ft(F) — ^1 - S. w .+S_.U). (4.9b) 

j~ij ^ 


. H , 

* ■ * . 




d.'. = d. . - -T d, , 6. . 
ij .ij 3 kk ij 


= 


1 


j “ 2 yax. ■ ax. 




fs S 
mn nm 

f 


0(F) = Plastic potential with 0(F) = 0 for F < 0 and 

T ~ — S S 
^2 ~ 2 mn nm 






To inclu(fe *the* elastic -vise opla Stic effects the constitutive equations, and 
conservation equations must be solved simultaneously with the aid of an equa- 


tion of state 


P = P(P,e ) 


(4.10) 
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4.2 FINITE ELEMENT ANALOGUE BY METHODS OF WEIGHTED RESIDUALS 

The finite element formulation posed herein is based on the conservation 
equations (4.5) through (4.7), the constitutive equations (4.9) and the Los Alamos 
equation of state. In principle, five conservation equations and nine constitutive 
equations must be solved simultaneously by the proposed finite element method. 
This approach, however, is not practical because it requires an extremely large 
storage space and possibly a large amount of computation time. Therefore, an 
alternative approach is introduced to solve each set of equations separately but 
with certain iteration procedures to couple the two systems of equations. ^ 

In this subsection, we discuss two finite element models constructed wil !i 
methods of weighted residuals, l.e., the Galerkin approach and the least squares 
approach. THe investigations were conducted to determine whether the finite 
element concept with conventional methods of weighted residuals can be applied 
to solve the high velocity impact problem, which is governed by a system of 
hyperbolic equations. Also they serve the purposes of program debugging and 
the development of a number of subprograms which are subsequently used in 
an improved finite element model. The developed procedures have been applied 
to solve a number of heat transfer problems, and impact problems as well. The 
findings from these numerical computations are presented and discussed in 
Section 6. 

4.2.1 The Galerkin Approach 

Consider a nonlinear boundary value problem of the form 

r (4.11) 

where a nonlinear (differential) matrix operator, A is the column vector 
of unknowns, and F is a column vector of the known quantities. Writing the 
conservation equations in this form, we have 
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(15 



^ + V. +0 

at ^ ^ ^i 

0 

0 


P 


& 

0 . 

s 

0 , » 
-KV“+ F- 
o t 1 

si; + ® - 0 

1 


V. 

J 


0 

0 

V. + 0 

9 t 1 o X. 

* ^ — 

« 

E 

L > 


and 


«: 0 * 


. s.. 


e. 


r = 


« 


I.SJ 


■ 

» 

I 


> 3 — (a. .) + F 

' ij j 


|-Q — (a. .V. ),+ F.v, 
i^3x. ' IJ j'i Eli 




(4.12). 


(4.13) 


Sim^larily, the conistitutive Eq. (4^t9a) can be written in the same form as 


. i « 




(4.14) 


and 


r.. = C.., , df- - S. ■ 03 • + S. 03 . . ' 

1 j 1 jk£ kfi im jm ixn im 


(4.15) 
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In the finite element/Galerkin method, approximations to p, V., and 


S. . are sought in the form 


P = N^Ps 


V = N V 

~ S ^8 


E = N E 

S 6 


s = 1,2, 


, n 


(4.16) 


S = N S 
s 


wherein n is the total number of unknowns, N are the shape (or trail)' func- 

s 

tions in space, p^, V^, E^, and^^ are the unknown time -dependent parameters 

to be determined; here V (tilda number) represents the vector V = jv^|, 
i = 1, 2, 3 and represents the deviatoric stress tensor S = J ^ 2, 3. 

(Note: The stress tensor is computed using the algebraic equation 4.8). in 
(4.16), we have assumed for simplicity and for computational convenience, the 
same type (linear, quadratic, or cubic) of approximation for all the variables. 


Upon substitution and applying the Galerkin technique, we obtain for the 
conservation equations 


K.| 

i'.l 



1 <P i 
1 1 St 

- 'o * 

- jO , 

(4.17) 

K.| 

|t| 

+ 

[®rs] 

1 (v 1 
1 r®i 

= (cl 

|~r| 

(4.18) 

[Arsj 

l-.l 

+ 

l^rs 

Ie I 

1 s) 

= [d { 
1 >•) 

(4.19) 


for r, s = 1, 2, . . . ,n, where the superposed dot indicates partial differentiation 
with respect to time, t, and 


A 


P3 


=4 N^N^dV. [0N^N^+v.N^ 

KrJ = | sl" “'ji + j 


dV 
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and 


D = 


/ 


N 


3x. 

1 


(a., V.) + F. V. } d¥ 

y 11 


for i,j = 1,2 and 3. The last two expressions in alternate forms can be re- 
written using Greenes theorem, 


= / <N^T.)dS+/ 


Bn 


a., + F. N 


ji j r 


dV 


D 


= / 


(N T.v.) dS + f (N F. - a.,|v. d¥ 

r J J ^ y r j Bx. yij j 


9N 


Equations (4.17) through (4.19) are the finite element analog of the conservation 
equations. '■ 


The finite element analog for the deviatoric stresses can be obtained in 
the similar way. Foir instance, witft Eq, (4.9b), we have 




(4.20) 


for r,s = 1,2, ...,n. Here 


i = /* N 
r s J r 


N d-V 

s 
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B 


rs 


■/ 


Zfi y N_ N_ + V, N 




r s 1 r 9x 


-^Id^, i = l,2,3, 


Q 


ki 


= /n^ 


> - Skm Vi + Vi “km k. t - = 1. 2. 3 


If the viscoplastic term is retained on the right hand side, then 


i = f V. 
rs J 1 


8N 


and 




= / 

jj- L 


- 70(F) 


kf 

yjhl 


- S w S . (J. > 

km mi mi km 


The conservation equations and constitutive equations thus formulated will 
be solved iteratively. The iterative procedures start with solution at previous 
time step and proceed as follows: 


1. Solve the conservation equations, Eqs. (4,17) through (4.19), to predict 
the field variables (p, 6 , u, v, and w) for next time step, 

2. Compute pressure using appropriate equation of state. 

3. Solve the constitutive equations (4.20) to predict the deviatoric stresses 
for the same time step. 

4. Use results obtained in Steps 2 and 3 and update the stress field. 

5. The above steps are to be repeated until a prescribed convergence 
criterion is met. If so, proceed to next time step. 
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4.2.2 The Least Squares Approach 

The possibility of utilizing the concept of least squares to both space 
and time has also been explored. To do this we write the equations to be 
solved in a symbolic form 

Q + A Q = L (4.21) 


in which Q represents the unknown function considered, L denotes the corre- 
sponding terms on the right-hand side, and A is defined as 


A = © 


+ ''k 


8 

9Xk 


(4.22) 


For a one-step scheme, the solution can be assumed in the following form 


■ 4 ' 


(n) 


Q - N, |(At -t) Q. + t Q 




/At 


(4.23) 


Again, N/s represent the space dependent shape functions, Q. and Q. 

^ (n ) (n 1 ^ 

are the unknown at nodal points for time equal to t' and , respectiv 


(rM 1) 
i vel y , 


with 0 < t < At. 


Upon substituting Eq. (4.23) into (4.21) with A ^ denoting A(NJ, the residual 
i s obtained as 


^ = A [- N, + (At - t) A J q/") + ^ (N. +tA.)Q.("+l) . L 

from which the weighting function can be obtained by evaluating the partial 
derivatives of R with respect to each Finally, the system of algebraic 

equations is obtained by setting 

ff 

tv ^ 
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from which a recurrence formula can be derived by carrying out the integration 
with respect to time explicitly. 

For instance, if and L are both assumed to be constant in the time 
interval and independent of the unknowns Then the weighting function is 


8R 




(n 


— ) 


The corresponding system of algebraic equations, with the time integration 
carried out, finally becomes 


where 




f I + A; N.) At (a. a, ) At^ 

K.j = J In. n. + '-i— J i-Jd — + v_L_iZ 




A. At 
1 


2 . /\^ 


(4.25) 


d^ 


or, equivalently, 


K.. AQ. = Li 


in which AQj is the increment of the unknown function at node j, 


AQ, = . Q (n) 

J 


j 


(4.25a) 


(2.24) 


4.3 ISOPARAMETRIC ELEMENTS AND NUMERICAL INTEGRATION 


As is seen, the use of finite elements discretizes a continuum problem 
and establishes a system of algebraic equations, whose coefficients are ex- 
pressed in terms of the products of element shape functions. The choice of 
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element type in finite element analysis is usually dictated by considerations 
of accuracy, computational efficiency, and the specific problem under study. 
Clough [28] has shown that for three-dimensional finite element analysis, the 
serendipity ’ elements with isoparametric formulation are superior to other 
solid element^ with respect to the above considerations especially in the 
present promem where an adequate representation of the geometry is essential. 
The term ’’isoparametric" means that the same shape functions are used to 
define both the geometry and the Crnknown function. In the present study, use 
will he made of the linear, quadratic, and cubic elements of the " serendipity" 

family. After performing certain transformation (mapping), these brick-type 

( 

elements will deform to yield curved surfaces. The following is a description 
of these elements. >. 

In the "serendipity" brick element, most of the nodes are located on ex- 
ternal edges. In fact, the linear, quadratic, and cubic elements contain no 
internal nodes at all. The corresponding shape functions are listed below 
using the notation of Zienkiewicz [29] . 

Let 

So = r, Ti. and 


where (^ , q , are the local coordinates and denote the coordi- 

nates of nodal points of the cube (see Fig. 4-1), 

Linear Element (8 nodes) 

The shape functions in this case are defined by 

- N. = |(i + e„)(i + ti^Hi + ^^), i = i,2 



8 
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Quadratic Element (2Q nodes) 

For corner nodes 

Nj -i(l + e„)(l + r,^)(l+ -2) 

Typical mid -side node 

g. = 0, ri- =+ 1, Cj = + 1 

N. = 1(1 - (1 + (1 + 

Cubic Element (32 nodes) 

For corner nodes 

- I"?] 

Typical side node 

^=+-3.r|.=+I, ^j=+l 

Ni = ^ (1 - e^)(i + 9i^)(i + ,i^)(i + 

Numerical Integration 

In finite element analysis, the matrices defining element properties, 

i 

e.g., stiffness, etc., must be found. These will be of the form 


I 


Iff 


[G(x,y, z)] dxdydz 


(4.26) 


in which the expression G depends on the equation being solved and the shape 
functions N. and/or their derivatives with respect to the global coordinate 
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system (x, y, z). As is seen earlier, the shape functions are written in the 
local coordinates T], therefore certain transformations must be per- 
formed so that Eq. (4.26) can be evaluated in the local coordinate system. 

First, the expression G(x, y, z), which^i, involves shape functions N. and 
■ their derivatives in the global coordinati^ sysi;em (x, y, z) must be transformed 
into those in the local coordinate system (|, When isoparametric formu- 

lation is used, the relationship between global coordinates and local coordinates 
is defined by 


X = N. (4,11,;) X. 

y = N. (4 ,t),;) y. (sum on i) (4.27) 

and 

z = N. z. 


in which are the shape functions as defined before and x^, z.^ represent 

the global nodal coordinates. By this transformation, the originally right 
prism in the (^, r|,C) space will become distorted in the (x, y, z) space. Now, 
since the shape functions are defined locally in finite element analysis, no 
transformation is necessary. However, the derivatives of shape functions 
with respect to (x, y, z) must be transformed into the local coordinate system. 
This can be done as follows. 

Using the chain rule, one has 
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where 




is the Jacobian matrix. Therefore, one obtains 


a N. 
1 

9 X 


9N. 

.Tf- 


9N. 
1 

9 z 


, I 


-1 


9N. 

1 


9N. 


9N.. 

1 


The*Jacobian matrii;^ |j| and its inverse, in turn, can be determined from 
(4,28). More explilgitly, the Jacobean matrix becomes 


1^] 




9Nj 

9N- 

A 

. . . . 

9Nj 




9 q 



9N • 
4*1. 


9 ^ 

a ^ 5- 



1 


^2 ^2 


(4.28) 


and inverse cari -be dete rmined,,.subse quentiy . 

Secondly, a transformation for the volume element must also be done, 


that is 


dxdydz =dc|f |jj d^drjd^ 


(4.29) 


® By combining (4.28) and (4,2^), one finally obtains, in place of (4.26), 

I « 

the fcillowing integral form ‘ | 

1 1 1 ■ ■ 

1 ‘I 

(4.30) 

L J / 

-1 -1 -1 


111 

■ ///• '[5(e. Ti, ^)j d| cirid^ 
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Thus the integration is carried out within the right prism and not in the 
complicated distorted shape. 

While the limits of the integration are simple in (4,30), unfortunately the 

___ 

explicit form of (G) is not. Therefore numerical integration usually has to be 

1 

resorted to. Essentially, (4.30) is approximated by the following form: 

^ ^m)| (4.30a) 

i = lj = lm=i L V J /J 

In the above, w., and are the weighting coefficients with G evaluatecj 

at the point (§ rj., C )• Herein, for simplicity, the number of integrating 
1 j m 

points in each direction was assumed to be the same. The numerical scheme 
presently used is the Gaussian quadrature, of which the points for evaluating 
G and hence the corresponding weights, are preselected to yield higher accu- 
racy with a fixed number of integration points (see Conte, pO]), Of course, 
within the isoparametric family, exact numerical Gaussian quadratures can 
be obtained if a sufficient number of points is used. 

The surface integrals can be expressed in a similar form 

I = JJ [H(x,y,z)] dS (4.31) 

S 

in which H(x, y, z) is obtained from the shape functions and/or their derivatives 
with respect to the global coordinate system (x,y,z) and S is a curved surface 
in space. Since the integrand is usually very complicated, expression (4,31) 
will be evaluated by numerical integration. 

First, the expression H (x, y, z) must be expressed in terms of the local 
coordinates analogous to the volume integrals, and to be evaluated 

on the appropriate surface. Second, a transformation for the area element 
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must be done so that the integration is performed on the surface defined by 
two of the local coordinates. For instance, the area element on a surface 
where t, is constant can be written as 


dS = det[j] dg drj 


(4.32) 


in which [jj is the modified Jacobian matrix evaluated on the surface con- 
sidered, The matrix, in turn, is defined as 


[j] = 


with 


3x 

h. 

6z 

H 

9? 

9e 

9x 

p. 

Bz 

9r? 

9r? 


CO ICO 
1 



w 


^ \9^) 


(4.33) 


The corresponding surface integral, througn the above transformations , 
thus becomes 

1 1 


■ ■ // n)] ^77 


(4.34) 


-1 -1 


Finally by using Gaussian quadrature, the surface integral is approximated 
by the following summation 


n n . 

^ = E S-i-j r)j)] (4.34a) 

i=l j=l 


in which H is evaluated at the Gaussian points (g., rj*), w. ,w. are the corre- 
spending weighting coefficients, and n denotes the total number of Gaussian 
points to be used in each direction. 
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4.4 TIME MARCHING SCHEMES 

As is seen previously, the Galerkin approach results in a system of 
ordinary differential equations with respect to time (see Eqs. (4.17) through 
(4.20)). These equations involve time derivatives of the primary conservative 
variables, p, V^, E and thus certain scheme must be chosen to integrate 
the system of equations in time to obtain the time history of the solution. 

As is seen, the equations being considered are all in the same form as 


a. . (f). + <l>. = y. 

J J-J J 


(4.35J 


One crucial point for the success of the present study is the choice of suitable 
time integration scheme to solve (4.35). There are numerous schemes avail- 
able; however, for the present problem, it is highly desirable to adopt some 
technique which is simple to apply, needs less storage locations, and should 
be numerically stable. We have investigated two such schemes as described 
in the following. 


Implicit Finite Difference Scheme; 
two consecutive time step solutions, 
(n+l)At. In particular, we assume 


In this scheme, use is made of the 
and , for time at nAt and 




(4.36) 


and 


+ (1 _ 0 ) 

J J J 


(4.37) 


0 < G < 1 
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Upon substituting (4,36) and (4.37) into (4.35) and rearranging, we 
finally obtain 

= + (l-e)p.j]^f> (4.38) 

Equation (4.38) is simple to apply and needs only one previous time 
step solution. Also, it can be shown for linear problems that the scheme is 
unconditionally stable for l/2 ^9^1 (see pi] ), thus one can avoid all 
stability worries and choose a value of At based solely on such considera- 
tions as desired accuracy. The above scheme, however, is only conditionally 
stable for 0 < 0 < l/2 and hence is not recommended for the problem pres- 
ently studied. 

Galerkin Process in Time; Alternatively, the Galerkin weighted 
residual process in time can also be applied to matrix equation (4.35) to 
obtain a recurrence relationship similar to (4.38). This approach allows 
a more comprehensive treatment and indeed possesses all the possible 
merits of the different variational processes suggested by Wilson and 
Nickell [32]. The recurrence could be written for several intervals simul- 
taneously thus necessitating more equations to be solved at each step but 
resulting in an improved accuracy and stability, thus allowing a larger 
time step to be used. 

Generally, within the interval we shall assume an interpolated form 
for each of the time dependent unlcnowns (f> ^ defined by its values at several 
time intervals 

0- = N^^\t)0^^^ (sum on i) (4.39) 

in which N^^^t) are appropriate shape functions defined continuously within 
the interval. 

For instance, if a linear interpolation is assumed, (4.39) can be 
written explicitly for the time interval (0 < t < At) as 



(n+ 1) 
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j • At ^ ^ At 


(4.40) 


From the above, the time derivative can be readily obtained as 


+ 




(n+1) 


:(4.41) 


i'n\ 

As the previous step solution (f)^ ' is known, only one weighted residual substi- 
tution needs to be evaluated. More specifically, Eqs.(4.40) and (4.41) are first 

e» ^ 

substituted into (4.35), then the resulting residual is multiplied by t/At, and the 
weighted residual is integrated over the time interval (0 < t < At). We obtain 
finally 








(4.42) 


It is seen that Eqs. (4.38) and (4.42) are similar and they become identical 
when 0 is set equal to 2/3 and is independent of time. The two approaches, 
however, are based on entirely different concepts. In the finite difference 
scheme time is discretized and Eq. (4.35) is solved directly; on the other hand, 
in the Galerkin approach the nodal unknowns are assumed to be continuosu func- 
tions of time and Eq. (4.35) is solved by some average process in time. The 
latter approach is therefore more general and can be conveniently extended to 
multiple time step procedures and higher order approximations in time. The; 
schemes discussed above represent only the simplest recurrence relations 
available and, of course, schemes of higher order approximation can be readily 
incorporated if circumstances require. 
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5* IMPROVED FINITE ELEMENT CODE BASED 
ON EULERIAN DESCRIPTION 

When shocks, compression waves or some other type of discontinuities 
occur in the material, the gradients of the dependent variables in the governing 
equations become very large in the neighborhoods of discontinuities. These 
large gradients lead the effective diffusion terms in the equations to be nega- 
tive in some regions, and cause numerical instabilities in the interaction pro- 
cess during the numerical computations, if conventional methods described in 
the previous section are used. To overcome these difficulties, an alternative 
finite element formulation was investigated. 

The present formulation consists of two primary portions. Firstly, the 
finite element formulation is constructed based on the theorem of weak solu- 
tions, so that the jump conditions can be satisfied automatically to take care 
of shock propagations. Secondly, a generalized two-step, time -splitting, shock 
smearing scheme has been developed and implemented. The scheme so con- 
structed has a capability to remedy the spurious oscillations arising due to 
numerical instabilities. 

In subsection 5.1, a general discussion on the theorem of weak solutions 
will be presented without proof. The governing equations in conservation form 
and the finite element analog of these equations, are presented in the sub- 
sequent subsection. The fomrulation of a two-step, time -splitting scheme is 
discussed in subsection 5,3. In order to predict the crater size and a numerical 
solution with adequate accuracy, the free surface of the projectile and target 
system must be treated with care. Discussion on this aspect is given in the 
last subsection. 


5-1 


LMSC-HREC TR D390900 


5.1 ON THE THEOREM OF WEAK SOLUTIONS 

The results from the theorem of weak solutions have been frequently 
applied in developing finite difference schemes for solving first order, non- 
linear hyperbolic equations. In most of these schemes, the advantages of the 
resulting forms deduced from the theorm are not fully utilized, since they are 
always in the integral form. This fa^t, ’^however , can be implemented in a 
simple way into a finite element scheme, and the advantages can be fully 
utilized. 


Various works on finite difference have indicated that, when the results 
of the theorem of weak solutions are to be used in the numerical scheme, the 
governing equations need to be converted into a conservation form (cf. Ritchrnyer 
and Morton [31] ). In what follows, therefore, we restrict ourselves in consider- 
ing the system of hyperbolic equations in conservation form. It is obvious that 
all the conservation equations can be written in this form. The constitutive 
equations, however, are not in the conservation form but can be recast into 
such form as shown in subsection 5.2. 


5.1.1 Theorem of Weak Solutions for First Order Quasilinear Equations 

Consider a Cauchy problem of a system of quasilinear hyperbolic equa- 
tions 


at 


G, k = 1, 2, 3 


( 5 . 1 ) 


k k 

on a cylinder R = Q (x) x T (t), where x = (x^, x^, x^), F = F (x, t, (j) ), and G - 
G (x, t, <^>). A class of piecewise smooth and piecewise continuous, vector val- 
ued functions (f) in R for t > 0 are called the weak solutions of (5.1) if the fol- 
lowing relation is satisfied. 

// ^ fj (5.2) 

R ' ^ ^ R 
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where ^ € C°® with compact support on R, i.e., ^ is any continuously 

differentiable function that vanishd*s on 9R, the boundary of R. ^ 

Lax [ 33 ], Oleinik [34], et al,, have shown tKht, if the following conditions 
are satisfied; 








1. For each k = 1,2,3,F has a continuQjus partial derivative 
with respect to (jt , and is bounded for (x, t) CR; 

m f I 


2. For bou4ded 


a2 pk 


and 


.2^k 


3(K 


34 , 


are dbntinuous. 


•c; I 


> 0 ; 


fill. 


0>i 

H* ■ 


3. The function G(x, t, (f>) has a continuous partial derivative 
with respect to <}>, ^ 


The generalized solution of Eq. (5, with a piece>^i,se continuous initial condition 
is unique, and satisMes Eq. (5,2), 

The proof of the theorem can.be found in the aforementioned references, 
and will not be presented here. HdPwever, we shall discuss in what follows 
some jimmediate coiisequences of the theorem directly related to our proposed 
formiSlation. 0 

fi ■ * .5 

. t 

V 

5.1.2 Jump Condition 

Following the Green's theorem, Eq. (5,2) can be deduced to a form 


//(♦ft 


R 


+ dt 
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-fU [C (x,t) (/i(x,t)] dS7 dt + // 

R dQ xT 

// 


dS dt 


(IQ dt 


(5.3) 


R 


where is the k component of the outward unit normal vector on 9 , the 
boundary of Q, 


Observe now the integrals 


I 


V 


aF df^dt 

R 


(5.4) 


for some t^cT, xef^; and 


1 

s 



(x, t) (x, t) 


dS dt 


dQ xT 


(5.5) 


for all xeOr^, and te T, Since ^ (x, t) is piecewise continuous function with 

compact support in R, it is obvious that I = 0 if there is no discontinuity in 

® 1 
R. Otherwise, we may deduce without difficulty the jump of the function F^ 

across the discontinuity, i.e.. 


, 3 Q xT 

s 

where is the surface of the shock layer, [[i^ denotes the jump of F^, 

the k^^ component of F, across the shock. 
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I 


For the integral if there is no discontinuity on R, Eq. (5.4) can be i;?. 


reduced to 


^ (x,t) </> (x,t) dfi - J“ ^ (x, o) «/> (x, o) df2 


Q 


Q 


Across a shock, however, the Leibnitz's rule yields 

^ I 

: M‘j 






(Xjsl) y c^(ti,|) dS(n) de 


+ 

T •, 




-y,4(X2.e) •/'(x^p?) y* c^(^.^) ds(ri) dg 


T 9 


9 an^(t) 


■,- « j 


• ^ 

wheij^ the subscripts 1, 2 denote the values at the upstream and downstream 

of th-fe shock, respectively. Here, t*he speed of sound is defined as 

9 '* 

I 




c « = 

s 


as^ 

9t" “k 


The last two terms in the above integral is no more than the jump of i.e.. 



(x, t) (x,t) d^2 - 


/ 


?>(x,o 


(x, o) dQ 


Q 


S 


# 

8 
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+ y* t, (x.t) Cg (x.t) f(^]] dsdt 


9r2 xT 
s 


Since the entropy condition requires that 


(5.7) 


% = 0 (5.8) 

The substitution of (5.6) and (5.7) into (5.3) then yields 


R 


= / 


t (x, o) 4> (x, o) d$7 - 



i G dS^dt 


(5.9) 


This consequence implies that (5,9) satisfies the jump condition (5.8) auto- 
matically. 


5.2 FINITE ELEMENT ANALOGUE OF WEAK SOLUTIONS 
Recall that we are considering a system of equations 


^ 8F^ ^ 

Bt 9x, 
k 


(5.1) 


k k 

where F = F (x,t, <p), k = 1,2,3 and G = G (x, t, ^). In this section, we shall 
formulate the finite element analog of equations in the form of (5.1). To do 
this, we have to write the conservation equations and the constitutive equations 
for impact problems in the above form- 
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5,2.1 Governing Equations in Conservation Form 

It is obvious that the conservation equatiohs (4.1) throughiX4.3) are in 
the conservation form. The constitutive equations (4.9), however, are not in 
the form of (5,1). Therefore, we recast the constitutive equations as 


as.. 3(v, S..) ^ / S.. , 

M. + - IK ^ 0S.. + Zix (d!. - y4. -H- 


at 




ylhj 


“ S. (jJ S . CJ. 

im m. mj. im 
J 


I 


where 0 = div v. *Here, an additibnal term 0S. . appears on th& ®right-hand 
side, which is considered as the forcing term. 

Thus, for high velocity impact problems, the governing equations can 
be written in the ccSiservation form (5.1) with 


ft 


4 -i 


9 


V.;j=;,l,2,3 


S^j! i, j = 1,2,3 




'^k®ij’ " 1.2.3 


; k,£ = 1,2, 3 




ti' 

ft.'- 
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r 


G 


V. 


0 1 


F.; j- 1,2,3 

F- V.; j= 1, 2, 3 ? 

J J 

[ ( 0 - 2liy<i> /V^) S. . + 2M d. . 


- S. 00 . 

im mj 




(5.10) 


5,2.2 Finite Element Formulation 

Suppose that the discrete approximation is constructed by appropriate 
interpolation functions in each element, and take ^ (x, t) defined in (5.9) to be 
the weighting function with compact support in the cylinder AR = x T. Here, 
is the volume of an element e, T = [ 0, At] and At is the time step in the 
time-split scheme. In order to utilize integral relation (5.9) to minimize the 
error, we have to approximate the (weighting) function ^ by the shape function 
at the previous time step, i.e., 




r 


(1 



1, 2, . . . m 


( 5 . 11 ) 


where m is the number of nodes in the element e. Also, approximate the 
solution to (5.1) as 




(n) 


JL 

At 


an+1) 


, s = 1 


m 


( 5 . 12 ) 
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With these expressions, we obtain a matrix equation in the form of 

(KJ - ^ = ([^rs] + ^ jc^j (5.13) 

where 

f ^^r 

B = / V, -r — ^ N d V 
rs J k 0x, s 

=/!«'' Sr (5.14) 

Here, the vector H = F - <^v. 

5.2.3 Remarks 

Instead of using Eq- (5.11), we may as well as approximate the wcightin}^ 

function by ^ N , r = 1,2,... ,m. In this case, however, tile integral 

r Zit r ^ 

formula (5.3) instead of (5.9) has to be used. Mathematically, these two form- 
ulations could achieve the same accuracy and possess a similar stability char- 
acter, Nevertheless, the surface integrals appearing in (5.3) would affect the 
numerical computations, as will be shown in the numerical experiments. 


This argument also explains the difficulties encountered in the conven- 
tional Galerkin method. It is essential that, after integration by parts, the 
Galerkin's formulation as shown in (4.42) is similar to the matrix equation 
(5.13). The only difference between (4.42) and (5.13) is their relaxation factors, 
i.e., a coefficient 2/3 in (4.42) is replaced in (5.13) by a factor -1/3, However, 
the forcing vector |^j.| (4*42) contains surface integrals as in the afore- 

mentioned case when C = Trr N and (5.3) are used. Moreover, matrix 

^r At r ' 
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equation (4.42) does not satisfy the jump conditions. Therefore, the con- 
ventional Galerkin formulation is considered less significant for solving 
hyperbolic equations. 

It is important to note that, however, neither the use of t = — N nor 
/ t \ r At r . 

= [l - ^1 can always represent accurately the class of test functions 
having compact support on R if the conventional shape functions are used. 

These known classes of shape functions may not vanish or equal to a same 
constant on all boundaries of the cylinder R, although we have presumed 
that they do as required by the theorem of weak solutions. The error arises . 
from this aspect, nevertheless, can be reduced easily by treating the boundary 
conditions with care during the numerical computations. Meanwhile, constructing 
a new class of shape functions does not seem to be practical due to the difficulties 
involved. 

5.3 GENERALIZED TWO-STEP, TIME-SPLITTING SCHEME 

It is well known that the finite element method with conventional assembly 
techniques is equivalent to the centered finite difference scheme. When the 
method is applied to solve flow problems containing shock waves, numerical 
instabilities generally arise primarily due to the lack of dissipative terms. 

In the finite difference approach, this difficulty can be overcome by either one, 
or a combination of the following schemes: 

• Introducing artificial viscosity, 

• Replacing the center difference by a noncocentric 
difference scheme, 

• Utilizing the two-step finite -difference scheme, or 

• Introducing Lax-Wendroff* s second order correction, 
etc . 

All of the above variations except the second have been adopted in the curreni 
finite element codes to solve the impact problem. It is found that the last tw > 
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approaches yielded promising results. However, due to the excessive compu- 
tation time needed with the Lax-Wendroff second order scheme, we have 
decided to adopt the two-step, time -splitting procedure. 

Consider a non-linear matrix equation 


([A„] - = ([V.l 


+ At 




where [A ], [B ] are square matrices, ic > is a column matrix, and ^ 
rs rs-' 

denotes the solutiom of at the n time step. In the present case, the factor 
6 = 1/3, the matrices 


(5.14) 

(n) 


i ; 
■I... 


n 


r 

^rs -/ dx 


n 


V dn 
, - k ■ . j 


r, s = 1.2 i 


and the vector 








H' 




where - Vj^<^ . 




Following the general approach of two-step procedure (see Richtmyer 
and Morton [31]), the solution of (5.14) can be solved by the following steps: 


(KsJ-a At = ([A^J + aAt(l-0)[B<"^>]) + aAt j (5.15) 
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(KJ-Ate = ([A^J + At(I-0) (5.16) 

where 

>+l) ^ aB<"; + (l-a) B|"/^> 

= aC<'’) + (1-a) 

Here a is a relaxation factor, and the constant a may be an integer or a frac- 
tion multiplier of the time step. As will be seen in Section 6, the numerical 
computations become more stable as the parameter a increases. The mathe-* 
matical aspect of the present formulation is yet to be studied in depth. In 
particular, the relation between a and the r^tio of the time step size and the 
element size must be determined in order to obtain an optimal accuracy for 
the general impact problems. 

5.4 FREE SURFACE CONSIDERATIONS 

In order to predict the crater size and produce a numerical solution 
with adequate accuracy, the free surface of the projectile and target system 
must be treated with care. Essentially, a free surface has zero pressure, 
but with its geometric shape changing with time. In the present formulation, 
the imposition of zero pressure on the free surface could be carried out con- 
veniently by proper consideration of the boundary integrals in the resulting 
algebraic equations; however, the adjustment of the free surface location is 
somewhat more complicated. There are several ways of handling the free 
surface problem. In this subsection, two approaches which are believed to 
be suitable for the present technique will be discussed. One of the approaches 
is to consider the equation governing the free surface motion as a part of the 
governing equations. The other one is simply to compute the solution on the 
free surface by an extrapolation procedure. 
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5.4.1 Method Accommodating the Free Surface Equation 


Consider the free surface in Fig. 5-1, which separates the interior of 
the target from the void region, and let h(x, y,t) be the vertical distance from 
the x-y plane to the free surface. An equation governing the free surface 
movement can be derived as follows. 

In Fig. 5-2, let A and B be two adjacent points on the free surface at 
time, t, with vertical distance defined by h(x, y, t) and h (x + uAt, y + vAt, t) 
respectively. After At, point A moves to point A' whose vertical distance is 
h(x + uAt, y + vAt, t + At). From the figure, we have 

Ah = h(x + uAt, y + vAt, t + At) - h (x + uAt, y + vAt, t) 

p = h (x + uAt, y + vAt, t + At) - h (x, y, t) + h (x, y, t) - h (x + uAt, y + vAt, t) 

= wAt - h(x + uAt, y + vAt, t) + h(x, y,t) 


Dividing through by At and letting At — ► 0, we have 

I 


dh 

8t 


V - u 


ah 

Sx 


v 


ay 


or equivalently, with first order approximation 


Ah 




(w - u 


3x 



(5.17) 


(5.18) 


As mentioned above, to analyze the free surface motion, we consider 
Eq. (5.17) as a part of the governing equations, and is computed successively 
with the conservation equations, the constitutive equations, and the equation 
of state. Equation (5.17Vjmust be solved to find h, or equivalently Ah, at 
each node. Here we confine our attention only to the vertical movement of 
the free surface. 
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Fig, 5-1 - Target-Projectile Configuration with Target Free 
Surface Defined by h(x, y, t) 


A' 



Fig? 5-2 - Motion of Points on Free Surface 
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The advantage of this approach is its formulation can be incorporated 
into the coupled Eulerian-Lagrangian code directly as part of the code in the 
rezoning process. However, since Eq.(5.l7) has the similar form as the 
governing equations, similar problems on numerical instability may arise 
during the computations. This problem remains to be investigated in future 
studies. 

5.4.2 Method of Extrapolation 

Recently, an improved numerical finite element procedure is proposed 
by France [35], The approach is applied to deal with two-dimensional steady 
state free surface problems. As an alternative to adjust the mesh in order 
to accommodate the movement and subsequent location of the free surface as 
in the previous method, extrapolation is used in conjunction with a fixed finite- 
element grid. The exact position of the free surface is determined when the 
imposed boundary conditions are satisfied. 

The extension of this approach into three-dimensional problems is 
straight-forward. The following illustrates the procedure with a linear brick 
type element. A typical element in the physical domain with the interface 
passing through is depicted in Fig. 5-3. Here the shaded area represents the 
free surface which cuts the edges of the element at a,b,c,d. The nodal points 
of the element are denoted by 1, 2, . . , 8. 

Let 4* t) ^ function whose value on the free surface is specified. 

At each time step, we assume that the function can be approximated as 

^ ~ ^-^>2,..,m (5.19) 

Then, the value of 4^ at the point a can be expressed as 

4^^= i [(I +h^)4-j Ki -h^)4-5J (5.20) 
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where (1, -I,h ) is the location of a based on the local coordinate system. 

SL , 

Thus, the nodal value of 4/ at node 5 can be obtained as 


= [2.];^ - (1 - h^)4-j]/(l +h^) (5.21) 

The values of 4* a-t nodes 6, 7 and 8 can be obtained in the similar way. The ' 
matrix equation for the entire boundary surface is thus generated in terms of 
the values at points of the free surface cutting the edges of all ” boundary" ele- 
ments. But all values of 4< at "boundary" nodes are computed from the previous;:^ 
time step, and v|j on a,b,c, . . . are all specified. Hence, the shape of the free 
surface can be obtained by solving this matrix equation. 

It is seen that this approach is simpler to formulate compared to the 
previous one, and the computations involve only algebraic equations. Since 
the element geometry is "fixed", it is only necessary to formulate the governing 
matrix equation for the initial cycle. These facts can naturally save cornputation 
time considerably. 
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5,5 SCHEME FOR SOLVING A LARGE SYSTEM OF EQUATIONS 

As is seen now, the numerical solution of the subject problem must deal 
with a system of governing partial differential equations. With the methods 
proposed herein, these equations are solved separately but coupled through 
iterations, thus reducing computer storage requirements to a great extent. 
However, for a general three-dimensional problem even with only one equation 
to solve, a solution with adequate accuracy will generally involve a large number 
of unknown parameters and hence result in a large system of algebraic equations. 
Therefore an effective scheme for solving such system of equations is obviously 
needed. 

The various methods to solve a large system of algebraic equations can 
be classified either as a direct elimination process or as an indirect iterative 
process. The direct elimination process solves the system at a minimum num- 
ber of arithmetic operations while storage is at a maximum. For a large scale 
problem, even with the banded nature of the matrix taken into consideration, the 
storage can still easily '^exceed the core memory available on most existing 
computers. Thus, any large system equation solver would have to make pro- 
vision for efficient transfer of data between core memory and auxiliary memory. 
The iterative process, on the other hand, requires only a minimum storage while 
the number of arithmetic operations, due to the iterative nature, is not definite. 
For matrix with diagonal dominance, the convergence is fast, but for any other 
matrix the convergence is usually slow or even not convergent at all. 

Traditionally, the finite element workers have favored the direct elimi- 
nation process for the following reason. The structural problems are mostly 
linear, or in case of nonlinear problems the nonlinear terms can conveniently 
be moved to the right-hand side with the load vector. Thus, for a given struc- 
ture, the decomposed coefficient matrix is invariant and the most time con- 
suming decomposition process needs only be performed once and for all. This 
decomposed matrix can then be used over and over to obtain a new solution 
corresponding to a new loading at a later time. It can also be used to obtain 
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an updated solution in an iterative scheme for nonlinear problems. This 
resolving capability of the direct elimination process makes it extremely 
attractive. 

Efforts have also been concentrated on storing only the nonzero coeffi- 
cients together with an auxiliary pointer matrix to record their locations. 

For input purposes, this scheme certainly saves a lot of storage. Unfortunately, 
nonzero coefficients are generated within the band in the decomposition process. 
Thus, the decomposed matrix would be relatively dense and it appears that there 
can be no saving on storage at output. In addition, some elaborate scheme is 
required to keep track of these generated nonzero coefficients. 

At this time, no definite conclusion has yet been reached whether to use 
the direct elimination process or the iterative process. It appears that storing 
the entire band matrix in core memory or auxiliary memory, and using the 
direct elimination process would be a good approach. Recently, two computer 
programs based on such an approach were published. The program due to 
Wilson et al. [36] was written for a positive definite symmetrical band matrix, 
while the program by Vendhan et al. [ 37] can take either symmetric or un- 
symmetrical band matrix. On the other hand, the iterative process such as 
the ^'Frontal Solution Technique” [38] is also very attractive for its minimal 
storage requirements and consideration of the nonlinearity of a problem. All 
these equation solvers are under study for possible inclusion in the final com- 
puter program. 
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6. TEST PROBLEMS AND NUMERICAL RESULTS 

In this section we discuss, in detail, both the inviscid hydrodynamic 
model and hydroelasto-viscoplastic model of a one -dimensional impact problem. 
Subroutines for numerical evaluation of shape functions, generation of matrix 
equations, as well as the equation solver, etc., are tested and debugged by solv- 
ing various cases of heat conduction problems. The finite element/Oal erkin 
procedure discussed in Section 4 is employed to solve the problem, and 
the details are presented in the Subsection 6*1. In the second subsection, 
we shall discuss the inviscid hydrodynamic code in some extent. Included in 
the discussions will be the numerical experiments of the given test problems 
in impacts, the comparisons between the methods of weighted residuals and 
the improved scheme discussed in Section 5, and some related aspects con- 
cerning the numerical technique. The hydroelasto-viscoplastic model will be 
discussed in the last subsection. 


6.1 HEAT CONDUCTION IN SOLIDS 

The heat conduction problem was chosen to test and debug a number of 
subroutines such as the numerical evaluation of shape functions and their deriv- 
atives for various three-dimensional isoparametric elements, assembly routine, 
time marching schemes, and equation solver, etc. The problem had been analyzed 
rather thoroughly to give us confidence in these subroutines, however, for brevity, 
only major results are presented herein. 

As is well known, the heat conduction problem is governed by the diffusion 
equation in the following form 


^ 

8t 


T 


( 6 . 1 ) 
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subjected to an initial condition and to a set of boundary conditions that are 
admissible to the equation. Note that the variables in (6.1) and in what follows 
are all normalized by the characteristic length, time and temperature. 

The steady state solution of the problem, i.e., the asymptotic solution 
of (6.1) at t — ► 00 , is identical to the solution of Laplace equation: 


V^T = ^ + ^ + 0 

9x^ 8y^ az^ 


( 6 . 2 ) 


subject to the same set of boundary conditions of either specifying the temper- 
ature, or the normal derivative of the temperature, or a combination of both, 
as in the unsteady case. This problem, when cast in a variational form, has 
the following integral expression 

- yy'(qT + iaT2)dS 

s 


It can be shown that, upon minimizing the above integral (i.e., 51 = 0), one will 
obtain (6.2) together with the following natural boundary condition. 


ar 

an 


q 1 aT 


(6.4) 
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Standard finite element procedures are then applied, based on (6-3) to 
obtain a system of algebraic equations in the form; 


|3.. T, = V. 
ij J 


where 


and 


6.. = /YT /"n. N. + N. N. . + N. N. \ d5,f 
V 

-jfy'aN. N. dS 


Vi = ff qNj dS 

s 


(6.5) 


( 6 . 6 ) 


In the above, (3^^ is the influence coefficient matrix, is the load matrix, 
and the repeated index j implies summation from 1 to m, the total number 
of unknown parameters. N. and N. etc., are the shape functions and their 

1 X , X 

derivatives, respectively. Since isoparametric elements are used, numerical 
integration scheme described previously is used to obtain these matrices. 


Two sample problems were chosen to check out the computer codes. 
They are the heat conduction in a cube and that in a hollow sphere, for which 
analytic solutions are available for comparison purposes (see Carslaw and 
Jaeger. [39], pp. 177-179, and pp. 230-231). 

6.1.1 Steady State Heat Conduction in a Cube 

Let us consider a cube (see Fig. 6- la) defined by 


0<x<a, 0<y<b, 
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subject to the following boundary conditions 

T = T j on X = 0 
T = T 2 on X = a 
T = 0 on other surfaces 


The analytic solution in terms of infinite series is 

T =14 ^ ^ [tj sinh f(a-x)+ T 2 sinh fx] sin sin 

TT m=0 n=0 (2m + I) (2n + 1) sinh la. 


in which 


|j -(2m+l)?r J ^ ^ p2n+^l)^ j ^ 


1/2 


For numerical computations, the element mesh shown in Fig. 6 -la is used 
together with a=b = c = l, T^=0 and = 1* To impose the boundary condi- 
tions on the surfaces of the cube, the nodal parameters on these surfaces are 
set equal to 0.0 or 1.0 accordingly. However, because the temperature is 
double -valued along the four edges of the surface x = a, some average scheme 
has been adapted to take care of this singular behavior, with the four corner - 
node parameters set equal to l/3 and the remaining nodal parameters on the 
edges set equal to l/Z. An alternate approach to resolve this difficulty is to 
use another grid with finer mesh arranged in the aforementioned region. With 
the present uniform mesh (two elements in each direction), computations have 
been carried out using both the cubic and quadratic isoparametric elements. 

Our results, as shown in Fig. 6-2, appear to compare well with the series solu- 
tion and the property of symmetry about the y'- and z'-axes was also observed. 
Both sets of finite element results show some oscillation about the series solu- 
tion» this waviness is believed due to the relatively coarse mesh used and the 
singular behavior of the boundary condition along the four edges on surface 
X = a. 
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b. Along y' - or z'- axis 

Fig. 6-2 - Connparison of Predicted Temperature Distribution in a Cube 
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6.1.2 Steady State Heat Conduction in a Hollow Sphere 

The hollow sphere is defined by a < r < b, subject to the ^boundary 
conditions 


T = Tj on r = a 


and 


T ^ *on r = b 


I 




The analytic solution is in the following form ^ 

a (b-r) T, + b (r-a) 

T = f ^ ^ 

r(b-a) 

Because the solution is a function of r only, a prism as* shown in 
Fig. 6- lb was taken' for computations. The inner and outer radii of the 
hollow sphere wer^ chosen to be 1.0 and 2.0, respectively. The boundary 
conditions imposedton the prism ^re * 

T = 0.5 on r = 1.0 

T - 1.0 on r = 2.0 


and 


d T 

= 0 Idh other surfaces 

a n 


Again, both quadratic and cubic i.s^oparametric eftments have b^on used in 
the computations, with results shown in Fig. 6-3. ■ As is seen, there exists 
excellent agreement between our predicted results and the analytic solution, 
even with only one element. Results obtained by using finer meshes are 
almost identical to the analytic solution. The high accuracy achieved by the 
present analysis is obviously att:jibutable to the ability of isoparametric 
elements to represent both the geometry and the solution functiJh accurately. @ 
Also, for the preserft problem, thejwell posed boundary conditions (no singu- 
lar behavior) make it easier to obtfiin an accurate numerical solution with 
relatively few elements. | 
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6,1.3 Transient Solution 


The time marching schemes discussed in Subsection 4.4 have been tested 
numerically for the three-dimensional transient heat conduction problem. 1 he 
governing differential equation now has the form 


3T 

8t 


ah a^T a^T 

2 t 7 

3y^ dz 


( 6 . 1 ) 


with initial and appropriate boundary conditions. By applying the Galerkin 
technique with respect to the space variables, a system of algebraic equa- 
tions, with an additional term involving time derivative, is obtained 


a. . T, + 6 . . T. 
iJ J ' iJ J 




( 6 .?) 


where 

V 




+ N. N. 

y y 


V 


and 





N. dS 
1 


s 



dV // aN.N.dS 
S 


( 6 . 8 ) 


The steady state solution of the hollow - sphe re problem was again selected 
as the testing case. The tests include different lime marching schemes, size 
of time step, various order of elements, and effects of initial conditions. Tlje 
solution is considered as reaching steady state when certain prescribed con- 
vergence criterion is satisfied. The one presently used is that for every 
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undetermined parameter the difference between two consecutive time steps 

.4 . 

must be less then 10 . All the testing cases except one (At = 2,0 in Fig. 6-4) 

satisfy this criterion in less than ten time steps. Figures 6-4 through 6-8 
show some of the results. 

Figure 6-4 shows the temperature history at the point r = 1.5 obtained 
by the implicit finite difference scheme with 0 = 0.5 (i.e., the Crank-Nichoison 
type), using two linear elements. Although the results show some oscillation 
about the steady state solution the scheme is obviously stable, regardless of 
time step size. The scheme of two-step Galerkin in time (equivalent to 0 = 

2/3 in the implicit finite difference scheme) was also tested for the same 
problem, with results shown in Fig, 6-5, Oscillation about the steady state 
is seen to be reduced significantly with improved convergence rate, thus 
demonstrating the merits of the Galerkin process in time. Two additional 
tests were conducted for cases with quadratic and cubic isoparametric ele- 
ments together with two-step Galerkin process in time. The trends of con- 
vergence for these cases are shown in Figs. 6-6 and 6-7, respectively. For 
this particular problem, no significant difference is observed in the paths of 
convergence with various numbers of elements. 

Results in Figs. 6-5 through 6-7 were obtained by assuming zero sold- 

|- 

tion throughout the entire field (including those on the boundary) at time t= 0. 

To see how the initial condition affects the solution path, non-zero boundary 
conditions were also applied from the onset (t = 0), and results are shown in 
Fig. 6-8. Some significant disturbance for the first steps is noticed, especially 
for cases with a large time step. The disturbance, nevertheless, was quickly 
damped out and convergent solution is still attainable. 

6.1.4 Radiation Boundary Conditions 

Subroutines for the boundary integral involving a curved surface has been 
written and was checked by analyzing the problem of heat conduction in a hollow- 
sphere subject to ’’radiation'' boundary conditions (see [34], p. 19). Instead of 
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Fig. 6-4 - Predicted Temperature History of a Typical Point (r = 1.5) 

of the Hollow Sphere with Various Time Step Sizes Jtwo linear 
elements. Crank-Nicholson in time with 6 = 0.5) 
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specifying the temperature on the inner and outer surfaces, the boundary 
conditions are now | 

= hj (Tj - T) at r = a 

g=h 2 (T 2 -T) at r=b 

The above expressions indicate; at r = a radiation from the medium at T^, 
and at r = b radiation into the medium at T^. h^ and h^ represent the ratio 
of surface conductance to thermal conductivity of the medium. The analytical 
solution for this problem is 

a^ [b^h^ - rlbh^ - 1)] + b^ [r(ahj + 1) - a^hj 

rfb^h^ (ahj + 1) - a^hj (bh 2 - 1)] 

For numerical computations, we choose 

a = 1, b = 2, hj = 2, h^ = 1, ^2 ” ^ 

Then the exact solution becomes 


T 6r - 4 
5r 

Figure 6-9 shows the predicted temperature distributions using only- 

two elements, linear, quadratic or cubic. The numerical results obtained 

by- using higher order elements and the analytic solution agree very closely; 

the results obtained by using linear elements, though less accurate, are 

generally good. Figure 6-10 shows the temperature history at the point r = 

1.5 with various element representations for the same problem, with Galerkin 

process in time. In all the cases the solutions converged to the steady state 

-4 

solution (convergence criterion 10 ) within 20 time steps. 
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S.2 IMPACT PROBLEMS WITH HYDRODYNAMIC MODEL 

j.2.1 Description of a Typical Impact Problem 

As pointed out in our earlier discussions, a hydrodynamic model is a 

>T;ood approximation in the early stages of th«^ high velocity impact process, 

luring which the pressure is comparable to the shear strengths of the target 

material. As a starting point we begin with a numerical solution to the hydro - 

lynamic equations with the inviscid adiabatic approximation. In this case the 

conservation equations, with a.. = -P6.-, become 

^ J 


at 


(p) = - 


ax. 


(pvj) 


at 


(pVj) = 


-^(P) - (pvjv ) + pf. 

I 1 


(6.9) 


^(pe) = - 


al“ - a!: (p«vj) +pf,v. 


view of (4.4), the above equations can be written in the alternate form 


+ V, 7^^ + Bp =0 


at i ax. 


^ (V.) + V. ~ (V.) + ^ + F. 

at J I 8x. J ) j 


( 6 . 10 ) 


w + "i air = a!r "i 
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Consider the impact o£ a cube on a simi-ihfinite body of possibley dis- 
similar material at a velocity of V^. For simplicity we can assume that the 
target is a semi-infinite cube of dimension four time larger than that of the 
projectile. Because of the symmetry about the z-axis, it is sufficient to 
analyze the quadrant bounded by positive x,y,z-axes (see Fig. 6-11). A 
typical finite element mesh of the project] le -target configuration is shown 
in Fig. 6- 12, which is generated by a program subroutine (MESH). The sub- 
routine generates nodal numbers and their coordinates in the target as well 
as in the projectile. 

Initial Conditions; V/e assume that the target is at rest at time t= 0. 

Let Sp denote the set of nodal points in the projectile and denote the set of 
nodal points in the target. Then S^- ^ ^t nodal points 

common to the target and projectile (i.e., on the interface). For t- 0^, we 
have the following initial conditions for density, ]>ressure, and internal energy. 





P: 

= ^op’ 

i c s 

p 


0 

i C S U S 

1 


P 

P. 

= -P , 

i C S 

1 

o’ 

o 


- 0, 

i C S US 

D 1 


where the subscript ”i" refers to Ihe i^^ node, and being the initial 

densities of the target and the projectile, respectively, and is the pressure 
at the interface, and should be calcuJated from the Rankine -Hugoniot relations 
for the target and projectile, 


(1 

^1 


Projectile 

= Pop 

niHL + ug = 

= Pi (1 - ni)/2p„p 


Target 

^2 ^ '’ot ^ ^2 

(1 - n^) = q (6-0.) 

C2 = P. (1 - n2)/2P„t 


6-20 





L.MSC-HREC TR D390900 



Fig. 6-12 - A Typical Finite Element Mesh of the Projectile-Target 
Configuration 
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Here the subscript "1” refers to the projectile and "2” refers to the target; 
q is the particle velocity at the interface and U is the shock velocity. We 
have assumed in writing (6.12) that the initial pressure and internal energies 
are zero. From (6.12) we obtain 


- q = 


[1 - r,i) 


(6.13) 


(1 - r,^) 


(6.14) 


The interface pressure is given when = P^. When the projectile and 

the target are of the same material the interface particle velocity, q, is 
given by 


q = VJ2 


(6.15) 


To calculate the interface pressure P^ at impact we proceed as follows: 
We have from (6.13) and (6.14)) 


ri2 = 1 - P< 


( 6 . 1 6 


For various values of we can calculate from (6.16), wherein the pres- 
sure Pj is calculated from the Hugoniot equation of state. 

Initial conditions on the velocity components u., v. and w^ in x, y and 
z -directions , respectively, are 


u. 

1 

= V, = 0, 

i 

C 

^t 


w. 

1 

= 0. 

i 

C 

s - s 

P o 

(6.17) 

w. 

1 

= -V 

i 

C 

s - s 

p o 


w. 

1 

= -q. 

i 

c 

s 

o 



where q is computed from (.5.14) (q is equal to V /2 for like metal impacts)^’ 


6-23 


LMSC-HREC TR D390900 


Another way is to use conditions at t=0", the time prior to impact, as 
initial conditions. To do this, every unknown is set equal to zero except 
density and those quantities related to impact velocity such as velocity, 
momentum and total energy in the projectile. These initial data may be dis- 
continuous, which are, of course, taken care of conveniently in the present 
formulation. This set of initial conditions has been used in most of the com- 
putations. 

Boundary Conditions; As mentioned earlier the pressure must be zero 
on the free surfaces. Since the target is semi -infinite, the material particles 
far away from the impact region are unaffected. Hence, the internal energy 
and velocities must be zero and the density must be undisturbed. These 
boundary conditions, together with conditions for the plane of symmetry, are 
sketched in Fig. 6-13. 

Remark on the Equation of State; As mentioned earlier, there are two 
forms of equation of state being frequently used in high velocity impact prob- 
lems. They are the Los Alamos equation of state and the Tillotson's equation 
of state, both obtained from experiment. Obviously values of computed pres- 
sure distribution in materials depend heavily on the particular equation of state 
being used, which in turn will affect the entire numerical solution. Some un- 
certainty apparently exists regarding the two equations as they do not appear 
to agree closely in general and each seems to have its own range of validity. 
For instance, the Los Alamos equation of state and Tillotson's equation of 
state differ substantially in the case of aluminum (see Figs. 6-14 and 6-15) 
for internal energy values below 7.5 Mb-cm /gm, and agree closely in the 

3 

range of 40 to 43 Mb-cm /gm, as seen in Fig. 6-16. Therefore, caution musi 
be exercised when a specific equation of stale is employed in the numerical 
computations, and it seerns that a more accurate equation of state is needed. 

6.2.2 One -Dimensional Impact Problem 

A one -dimensional impact problem is used herein to test the validity 
of the formulations presented in Sections 4 and 5. The problem under 
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Internal Energy, pe (Mb) 

Fig. 6«14 - Variation of Pressure with Internal Energy (Aluminum 
= 2.702 gm/cm3pt = 0 to 1.5 Mb) 
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Fig. 6-15 - Variation of Pressure with Internal Energy (Aluminum 
p„ = 2.702 gm/cm3,p« = 2.7 to 7.5 Mb) 
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consideration is that an infinite plate traveling at a specified speed, v^, hits 
normally on another infinite plate at rest, as depicted in Fig.6-17a. The two 
plates are 30 cm thick, and composed of the same aluminum niaterial with 
initial density = 2.77 gm/cm^. The impact velocity is chosen to be = 
0.008 cm/jUsec for which a solution in elastic range of impact is available for 
comparison [40]. 

The problem was solved by the three-dimensional models described in 
Sections 4 and 5. The mesh and nodal numbering are generated by subroutine 
MESH and depicted in Fig. 6- 17b. The Rankine -Hugoniot pressure, is 
readily obtained, using the Los Alamos equation of state, 

P = 5.99868 X 10 ^ megabar (6.J8) 

s 

The following numerical computations are performed based on inviscid assump 
tion together with the Los Alamos equation of state. 

Solution Computed by Methods of Weighted Residual; It is found from 
the numerical solution of (6.10) that the Galerkin approach does not perform 
satisfactorily. Although the scheme appears to be stable up to the time com- 
puted, there is a sign that the pressure development in the material tends 
to grow indefinitely. This phenomenon becomes more severe when the impact 
velocity is increased. Spurious oscillations with fairly large amplitudes always 
gather near the wave front, and the peak pressure is seen to exceed the Hugonio 
pressure by about 50% to 100% as the velocity v^ increases from 0.008 to 0.75 
cm/p-sec. As remedies to the instabilities, some modifications such as incor- 
porating the mid-point Runge-Kutta scheme into the code, as well as other 
methods of weighted residuals such as the least squares approach are intro- 
duced. Figure 6-18 shows the time history of pressure at interface computed 
by the method of least squares. The results computed by a modified Galerkin 
approach using mid -point Rungc-Kutta scheme are found to be almost identical 
to those shown in Fig. 6-18. There is some improvement on the development 
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history of the pressure field. However, the pressure at the interface still 
shows a tendency to diverge as time increases. The numerical instabilities 
can also be detected by the spurious oscillations near the shock fronts as 
shown in Fig. 6 - 19. 

Solution Computed by Weak Solutio n Formulation: The pressure history 
at the interface computed using (5.15) and (5.16) with a= 1/2 is depicted in 
Fig. 6-20, and the corresponding pressure distributions at various times are 
shown in Fig. 6-21. It is seen that the results are improving compared to the 
ones computed by methods of weighted residuals as depicted in Figs. 6- 18 and 
6-19. The solution may be improved even further when the parameter a in 
(5.15) is increased. Figures 6-22 and 6-23 show the pressure history at the 
interface using a=2, and 4, respectively. The corresponding pressure dis- 
tributions are plotted in Figs. 6-24 and 6-25. 

All results discussed so far and depicted in Figs. 6-18 through 6-25 are 
computed using 16 linear elements. As shown in Figs. 6-21, 6-24 and 6-25, 
the overshoot due to the numerical instability decreases rapidly, and the zig- 
zag behavior near the wave front becomes less severe as the parameter a 
increases. These facts indicate that the numerical dissipation introduced in 
the scheme is proportional to the parameter a. 

The plots also show some substantial differences in amplitudes of the 
pressure waves among different a and a values. In addition, as shown in 
Figs. 6-22 and 6-23, the numerical dissipation increases as a decreases, 
and the phase lag becomes more apparent at the same time. These imply 
that various errors may enter in the computations if the scheme is over- 
dissipated. The sensitivity of the solutions upon the factor a, however, can 
he made inert by refining the mesh as shown in Figs. 6-26 through 6-29 where 
30 even spaced linear elements were used in the computations. In this case, 
the same time step size as employed in the 16 -element case was used. This 
verifies the statement mentioned in Section 5.3, namely, the factor a is indeed 
related directly to the ratio of the time step size and the element size. As the 
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Fig. 6-20 - Pressure Development at the Intensity Interface with 16 Linear Elements (a = 0.5, (X - 0.0) 
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-—— a - 2.00, a = 0,0 
a = 2.00. a - 0.25 



Fig. 6-26 - Pressure Development at the Interface Using 30 Linear Elements 
(a - 2.0, a = 0.0 and 0.25) 
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Fig. 6-27 — Pressure Devrelopment at the Interface Using 
30 Linear Elements (a - 4.0, a = 0.0) 
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mesh size reduced, the effect of a on the accuracy of the solution becomes 
less important. The results using 16 quadratic elerrients also show the same 
trend. In Figs. 6-30 and 6-31, the pressure history at the interface is depicted 
using a =2.0 and a =4.0, respectively, with different Ct. As expected, these 
results are improved compared to the cases with 16 linear elements, but are 
less smooth when compared with that using 30 linear elements. Similar com- 
parisons can also be made for the pressure distributions (see, e.g.. Figs. 6-24, 
6-25 for 16 linear elements, Figs. 6-28, 6-29 for 30 linear elements, and Figs. 
6-32, 6-33 for 16 quadratic elements). 

As is seen, all the cases seem to underpredict the peak pressure com- 
pared to the analytic solution. However, with the same given conditions and 
Los Almos equation of state, the pressure computed by Rankine-Hugoniot 

relations is P = 5.99868 x 10” megabar. This value is approximately 
s 

equal to the average value as shown in Figs. 6-28 and 6-29 which indicates 

that our numerical results are quite reasonable. Accordingly, P should be 

s 

the upper bound for the average pressure distribution in the material under 
high velocity impact, since the physical as well as numerical dissipative 
effects may actually reduce the pressure buildup in the material. Therefore, 
it is believed that the deviation of the analytic solution from the computed 
Rankine-Hugoniot pressure and from our results could be due to one of the 
following: 


• The assumptions made in the theoretical analysis 
lead to a too simplified model, so that the analytic 
solution was overestimated, 

• There are possible misprints related to the given 
conditions in the report we obtained, and 

• The analytic pressure distribution, Ps = 7.22 x lO 
was normalized by the impact velocity, Vq= 0.008 
cm//lsec, but was not mentioned in the report. 


These possibilities hopefully can be cleared up as soon as we locate the 
original report. 
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The speed of sound in a material can be obtained in the relation 



s 


={^) 

'ep 's 


Following the Gibb’s equation, one has 


(6.9) 


^2 p ap ^ ap 


( 6 . 10 ) 


where 6 is the specific internal energy which is now being considered as a 
function of total specific energy and the particle velocity. Thus the propagation 
speed of the pressure waves can be computed by substituting the nodal solution 
of the conservation equations and equation of state on the clement containing th'- 
wave front into (6.10). The numerical values of in the target obtained using 
Los Alamos equation of state at some typical times are shown in Table 6-1. 

It is seen that the values oscillate around the shock speed computed by Rankinc'- 
Hugoniot relation, i.e., = .53.569 cm/ (isec . In elasto-dynamic theory, 

on the other hand, the souri^ speed is represented by 


C 

s 




(1 - t^)E 

(1 + i') (1 - 2l/)p 


( 6 . 11 ) 


where i; is the Poisson ratio, E is the Young's modulus. With U = 0.33, E = 

10^ psi for aluminum, Eq. (6.11) gives C .= 0.60729 cm/ /isec, which deviate s 

®Ei. 

from the values obtained by the present approach by about 10%, This implies 
that the compressibility effects play an important role in the dynamic response 
of materials under high velocity impact loads. In addition, although in the 
present hydrodynamic model the material is assumed to be inviscid, the dis- 
sipation resulting from the equation of state and numerical viscosity as well 
may affect the local speed of the pressure waves. 


6-49 


-50 


Table 6-1 


PROPAGATION VELOCITY OF PRESSURE WAVES IN TARGET 
AT VARIOUS TIMES AFTER IMPACT 
(WITH 30 LINEAR ELEMENTS) 

Poisson Ratio, v = 0.33; Young's Modulus, E = 10^ psi; 
Impact Velocity, V^ = - 0.008 cm/jLLsec 


Time, t(jLLsec) 

0 

5.0 

10.0 

18.0 

26.0 

31.6 

35.0 

Shock Speed, 

C (cm//jLsec) 
s 

.53442 

.54123 

.53846 

.53880 

.53679 

.53801 

.53782 


Rankine-Hugoniot — C = .53569 cm/ fi sec 

s 


3-D Elastic Waves — = y (1 -f u) (l ^ em/jisec 
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Momentum and Energy Distributions; Momentum and total energy dis- 
tributions at t=20 and 30/lsec are depicted in Figs. 6-34 and 6-35 for 16 linear 
elements. Figs. 6-36 and 6-37 for 30 linear elements, and Figs. 6-38, 6-39 for 
16 quadratic elements. In these plots, the parameter a= 2.0 are used. A 
better accuracy using 30 linear elements against 16 quadratic elements is 
obvious. In particular, as seen in the rarefaction region, the momentum and 
total energy distributions computed using 16 quadratic elements are severely 
distorted by the numerical instabilities. Hence, a larger parameter a is 
needed, although it is not so obvious from the plots of pressure distributions. 
Depicted in Figs. 6-40 through 6-42 are momentum and total energy distribu- 
tions using various types of elements with a= 4.0, and as expected, better 
results are obtained. 

6.2.4 Three -Dimension,al Impact Problem 

3 

A simple problem of a 4 cm aluminum (p^ = 2.702 gm/cm ) cube im- 
pacting at a velocity of 2,6 cm/p sec on a semi-infinite cube (16 cm cube) 
was also tested. The numbering of nodes in the mesh is depicted in Fig. 

6-43, with the corresponding boundary conditions given in Fig. 6-13. 

For this problem, only preliminary results by the Galerkin procedures 
are available, which are shown in Figs. 6-44 and 6-45. Figure 6-44 shows 
the pressure variation with time at various nodes on the interface, while 
Fig. 6-45 depicts the pressure variation with distance into the target at various 
time. In the computations, Los Alamos equation of state was used and nega- 
tive pressures at any node in the projectile -target configuration were not 
allowed. Nor was the movement of free surface accounted for at the time 
of computation. Like the one -dimensional problem, the Galerkin formulation 
again indicated some numerical instabilities which must be remedied. 

Computations using the weak solution formulation is in progress, which 
will be discussed in the final report. 
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6.3 IMPACT PROBLEM WITH HYD ROE LA ST O- VISCOPLASTIC MODEL 

The governing equations with their finite element analogues are given 
in Section 5, The computer code will, at the end, be capable of handling large 
deformations, anisotropic materials, plastic yielding knd material fracture. 

The one -dimensional impact problem discussed earlier is used also to 
test the present model. The problem description and the finite element mesh 
have been shown in Figs.6-17a and 6-17b. In the computations, the material 
constants are assumed to be: 

3 

Initial density, = 2.77 gm/cm 
Shear modulus, p= 0.276 megabar 

The impact velocity, = 0.008 cm/jLLsec, and Los Alamos equation of state 
was also used for the present problem. For the purpose of comparisons, 
the problem was solved by both the Galerkin method and the weak- solution 
fo rmulation. 

Results Computed by Galerkin* s Procedures: Figures 6-46 through 6-48 
show the pressure and axial stress histories at the interface, the normal stress, 
momentum and total energy distributions computed by the Galerkin method. 
When they are compared to the results computed from the inviscid hydrody- 
namic code, it is clear that the spurious oscillations behind the shock front 
are smaller. However, the pressure development at the interface seems again 
to grow indefinitely, though the results computed up to the time is still finite. 
The axial stress, on the other hand, is found to be much lower than the Hugoniot 
pressure, P . Moreover, the normal stresses on planes perpendicular to tht 
axis are mostly positive, which imply the materials are under tensile stress 
in the axial direction. This obviously is not physically correct. These phe- 
nomena, therefore, indicate that the numerical results computed by Galerkin'. s 
method are totally unacceptable. 
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Fig. 6-46 —Axial Stress and Pressure Developments at the 
Interface Using Galerkin's Method (16 Linear 
El ements) 
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Results Computed by Weak Solution Formulation; The numerical solu- 
tion of the problem is much improved when the technique based on weak 
solution formulation is used. In Fig, 6-49, the pressure as well as the axial 
stress history are depicted against time. The parameters used in this com- 
putation are a = 4.0, a = 0.0, and the time step. At = O.ZjLtsec. When these 
results are compared with those in the inviscid case (see Figs.6-2Z and 6-23) 
some deviation is observed. In the inviscid hydrodynamic model, the pressure 
history computed using both a= 2.0 and a= 4,0 is found to be damped out signif- 
icantly as time increases. Using the present method with a = 4.0, however, 
results seem to suggest that less dissipations are present. The reason for 
this is not clear at the moment. Further mathematical study of this algorithm 
is needed to clarify the question. 

The distributions of the normal stresses at t = 10.0, 20.0 and 34.2/isec 

are depicted in Figs. 6-50 through 6-52, respectively. As indicated in the 

plots, the compressive axial stress a has an average value approximately 

zz ■ 

equal to the Hugoniot pressure, P , after some time of the impact. This, 

s 

again, suggests our numerical results are quite reasonable. For this partic- 
ular example, the impact velocity is moderate, and the viscous effects may 
come only from the equation of state but not from the constitutive equations, 
i.e., the shear stresses computed from the constitutive equations are negligibly 
small compared to the normal stresses. This implies that, the dynamic re- 
sponse of the materials due to impact is in the elastic range, and the stress 
wave must have an average value approximately equal to P^. 

From Figs. 6-50 through 6-52, we may also find the fact that, even the 
parameter a used in the computations is set to be 4.0, the spurious oscillations 
behind the wave front are still apparent. The lack of numerical dissipative 
effects in the scheme can also be detected in the interacting region of rare- 
faction waves as depicted in Fig. 6-53. To smooth out the zig-zag behavior, 
an increase of a can be employed, so long as the stability criterion of the 
present time marching scheme is not violated. A better way for reducing 
these oscillations, however, would be to refine the mesh and/or use higher 
order elements. Further studies on this aspect arc underway and detailed 
discussions will be presented in the forthcoming final report. 
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Fig. 6-49 “Axial Stress and Pressure Developments at the 
Interface (16 Linear Elements; a - 4.0, 0.0) 
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Fig 6-59 — Nornial Stresses Distr ibutions at t = 10.0 /isec (16 Linear 
Llemcnts; a ^ 4.0, O' 0.0) 
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Fig. 6-51 -Normal Stresses Distributions at t = 20.0 ^sec (16 Linear Elements; 
a ^ 4.0. a -r 0.0) 




LMSC- HREC TR D390‘K)0 






l.MSc: * MKl'X: 'i‘K 


7. SUMMARY AND DISCUSSION 

The objective of the current study is to develop a finite element com- 
puter program for the numerical solution of three-dimensional high velocity 
impact problems, based on the Eiilerian hydroelasto-viscoplastic formulation. 
Two models, namely, the inviscid hydrodynamic model and the hydroelasto- 
viscoplastic model were formulated. The theoretical basis and detailed form- 
ulations for the subject problem were discussed in the preceding sections. 
Computer programs based on the methods of weighted residuals and the thoorm 
of weak solution have been coded and debugged by running a number of test 
problems , 

For general impact problems, the conventional methods of weighted 
residual were found to be unsatisfactory. The main reasons are: 

• The formulations, such as in the methods of Galerkin 
and least squares, do not generally satisfy the jump 
conditions. 

• The finite element analog of the problem always lacks 
of dissipative terms, and thus causes severe numerical 
instabilities behind the discontinuity such as shocks, as 
well as in the interacting region of rarefaction waves. 

These difficulties can be overcome by the two-step, time- splitting finite 
element formulation based on the theorem of weak solutions. The success 
with this formulation has been demonstrated by a number of numerical experi- 
ments. As indicated by the numerical results, though the relaxation factor a 
is related directly upon the ratio of the time step size to the mesh size, the 
sensitivity of the solution upon a is rapidly reduced by refining the mesh and 
by using higher order elements. The value of the parameter a in Eq. (5. 15) can 
also be reduced by similar treatments to acquire the same stability criterion. 
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It is believed that, an optimal choice for a would be the satisfaction of the 
following condition: 

where c^ is the local sound speed, is the mesh size along Xj^- direction. 

Despite the successfulness of this presently developed technique, its 
mathematical structure needs to be investigated further. At the present time, 
the suitable choices of a and a rely mainly on numerical experiments ; the 
stability criterion for the present scheme has not been analyzed. These prob 
lems can only be resolved by a thorough mathematical analysis of the present 
procedure. 

In view of the findings up-to-date further work will be continued as 
follows: 

• Investigate the mathematical structure of the two-step, 
time -splitting, weak- solution formulation developed in 
this study. 

• Incorporate a subroutine to properly account for the 
free surface movement. 

• Incorporate a large system equation solver. 

• Include subprograms for plastic yielding and material 
fracture. 

• Develop a code for coupling Eulerian and Lagrangian 
modes. 

• Prepare and run the demonstration problems. 
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